Meteor generated plasma columns in E-region ionosphere: fields and diffusion 
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A meteoroid penetrating the Earth's atmosphere leaves behind a trail of dense plasma in the E-region iono- 
sphere, a region where electrons are strongly magnetized while ions are demagnetized due to their frequent 
collisions with neutrals. While radar measurements of meteor trail evolution have been collected and used to 
infer meteor and atmospheric properties since the 1950s, no accurate quantitative model of trail fields and diffu- 
sion exists. This paper describes a theory and simulations of trail plasma physics which applies to the majority 
of small meteors. Unlike earlier papers, our theory assumes a significant angle between the geomagnetic field 
and the plasma trail and includes the important interaction between the trail and the background ionospheric 
plasma. This study provides quantitative knowledge of the spatial distribution and dynamics of the plasma den- 
sity and electric field. This should enable meteor and atmospheric researchers to more accurately interpret radar 
observations of specular and non-specular meteor echoes. 
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I. INTRODUCTION 

Meteoroids impacting the Earth's upper atmosphere have a 
number of important consequences by: depositing material, 
most notably metals and dust; damaging spacecraft; creating 
layers of material which radars and ionosondes detect; modi- 
fying the plasma density and conductivity of the lower iono- 
sphere; leaving plasma columns which can be used for meteor 
burst communication and to monitor atmospheric conditions 
in the lower thermosphere. While large meteoroids gener- 
ate spectacular optical displays, particles much smaller than a 
sand grain comprise the majority of all meteoroids and repre- 
sent the major source of all extra-terrestrial material deposited 
in the Earth's atmosphere |lL|2|]. Such small meteoroids can 
rarely be observed by a naked eye or even by a sensitive opti- 
cal technique, but radars can easily detect them. Such obser- 
vations become possible because meteoroids frequently enter 
the Earth's atmosphere with a sufficient speed (11-72 km/s) 
and energy to cause the formation of a dense plasma visible 
to radars. The front edge of this plasma can often be observed 
by high-power large-aperture (HPLA) radars, while smaller, 
lower-power radars can detect the residual trail as either spec- 
ular or non-specular echoes. 

Specular echoes, usually observed by small-aperture and 
HF- VHF radars, originate from parts of the plasma trail where 
the wavevector of the backscattered wave is nearly perpendic- 
ular to the axis of a slowly diffusing, quasi-cylindrical plasma 
trail (for review, see |2]). Non-specular meteor echoes ob- 
served by HPLA VHF or UHF radars like the one shown in 
Fig. ^ originate from trails where the specular condition is 
not necessarily satisfied but where the radar typically points 
close to perpendicular to the geomagnetic field O B S 01- 
Non-specular echoes appear to result from small-scale elec- 
tron density irregularities caused by plasma turbulence and 
measured when the radar wavevector lies parallel to the ir- 
regularity wavevector but with a wavenumber twice the irreg- 
ularity wavenumber, satisfying the Bragg condition. A polar- 
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FIG. 1 : Non-specular radar echo (Jicamarca Radio Observatory near 
the magnetic equator, July 12, 2005, 3:43 AM local time). 



ization electric field, resulting from the plasma trail ambipo- 
lar diffusion, drives field-aligned instabilities which generate 
these irregularities 

Modeling specular echoes requires knowledge of the spa- 
tial and temporal distribution of the plasma trail density, while 
modeling non-specular trails also requires knowledge of the 
evolution and structure of the polarization electric field which 
drives the instabilities. In a series of papers 

USES HUE, 

Oppenheim, Dyrud and others have simulated the develop- 
ment of plasma instabilities in the meteor trail under con- 
ditions when the axis of an axially-symmetric cylindrical 
plasma trail was perfectly aligned with the geomagnetic field 
[or, in a two-dimensional (2D) case, making the trail a slab 
instead of a cylinder]. This is, however, a degenerate case be- 
cause most of meteor trails are oriented at a significant angle 
with respect to the magnetic field. In the general case, elec- 
trons have anisotropic responses to the electric field, which 
significantly complicates theoretical description of the dense 
plasma trail evolution. 

A number of studies have looked at the evolution of a dense 
plasma column in a collisional magnetized plasma. The ear- 
liest simulations had restricted box sizes which 
could not properly describe the actual ionospheric situation. 
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Furthermore, numerical simulations alone cannot provide the 
parameter dependence of the plasma density and polarization 
electric field, which is of importance for modeling plasma in- 
stabilities. At the same time, analytical theory, even an ap- 
proximate one, could provide such dependencies and would 
be an indispensable tool for the accurate interpretation of radar 
observations, meteoric and ionospheric diagnostics. 

The analytical theory by Jones 11711 represents a significant 
step forward. He proposed a 2D self-similar solution (SSS) 
of the meteor trail evolution, starting from an initial line den- 
sity. Using a combination of a Gaussian spatial distribution of 
the plasma density with a parabolic distribution of the electric 
potential, Jones developed a mathematical scheme, which de- 
scribes the initial evolution and structure of the trail plasma 
density reasonably well. However, it improperly describes the 
spatial structure of the polarization electric field because the 
assumed plasma density structure requires that the field goes 
to infinity with distance from the trail axis. As we show in 
this paper, that theory also fails to predict significant devia- 
tions from the self-similar diffusion at a later stage of the me- 
teor trail evolution. The reason for this is that the SSS does 
not account for interaction with the background ionospheric 
plasma. 

In the earlier stage of plasma trail diffusion, the background 
plasma density is usually small compared to the plasma den- 
sity within the trail. However, this low-density plasma plays 
a crucial role for carrying electric currents originating from 
the trail during its ambipolar diffusion. This current signifi- 
cantly affects the evolution and structure of the trail density in 
the later stage. Further, the electrodynamic interaction of the 
meteor trail with the background ionosphere provides a nat- 
ural restriction for the polarization electric field which drives 
plasma instabilities. The research presented in this paper in- 
cludes the interaction of the trail with the background iono- 
sphere. 

Note that in a recent theoretical paper on meteor trail diffu- 
sion, Robson filil (see also O) attempted to revise the previ- 
ous theories. Robson's approach, however, seems invalid be- 
cause, in addition to the quasineutral assumption, Robson as- 
sumed that the electron and ion fluxes along the magnetic field 
are equal. In the general case, the latter assumption causes the 
resulting electric field to have a significant non-zero curl. Es- 
timates show, however, that the contribution of the induction 
electric field to the trail dynamics in the E region is negligible, 
so that to high accuracy the electric field within and around the 
trail should be electrostatic, i.e., curl-free (see Sect. lVIlII ). 

In order to improve modeling of non-specular trails and en- 
hance our knowledge of plasma column diffusion and fields, 
we revisit this problem. As in the majority of the previous 
studies, we restrict ourselves to plasma columns that remain 
homogeneous along their length and assume constant electron 
and ion temperatures. Homogeneity along the trail implies 
that diffusion occurs only in the plane perpendicular to the 
trail, though electron fluxes have all three vector components. 
We further assume that the trail axis is directed at a sufficient 
angle to the geomagnetic field (in this paper, we restrict the 
analytical treatment to the purely orthogonal case). The paper 
differs from Jones 11711 in that it includes the important effect 



of the background ionospheric plasma. Finally, we consider 
only the case when there is no external electric fields or strong 
neutral winds which often exist in real E-region ionosphere, 
especially at high latitudes and in the equatorial region. We 
will describe these important effects in future papers. In this 
paper, we present results of both numerical simulations and 
analytical theory. The main result of this work is a quantitative 
description of meteor trail evolution and the polarization elec- 
tric field associated with its ambipolar diffusion. Results of 
our analytical model and simulations agree remarkably well. 
We expect that applying our theoretical model to radar echoes 
from trails should help researchers obtain useful information 
about meteoroids and the surrounding atmosphere. 

The paper is organized as follows. In Section|ll| we discuss 
qualitatively the ambipolar diffusion of meteor trails in the 
E-region ionosphere, introduce the concept of the residual po- 
tential and discuss some restrictions. In Section|III| we present 
the equations for meteor fields and diffusion based on a sim- 
ple two-fluid model model of a highly collisional isothermal 
plasma and formulate proper boundary conditions. In Section 
IIVI we discuss results of our 2D finite-element simulations, 
which give us useful insights into finding the proper approach 
for the analytical treatment of the problem. In Section [V] we 
present our analytical theory for the specific case of mutu- 
ally orthogonal meteor trail axis and the geomagnetic field. In 
Section IvTl we summarize our analytical results. In Section 
IVIII we compare our theory with simulations quantitatively. 
In Section fVlIII we discuss some caveats and implications of 
our theory. In Section Hxl we give a summary of the paper. 
AppendiceslAlto IFI give mathematical details of our analytical 
theory. 



II. QUALITATIVE DESCRIPTION OF METEOR TRAIL 
DIFFUSION 

In this section, we discuss qualitative aspects of the meteor 
trail diffusion which will allow us to formulate the problem 
and make appropriate simplifications. 



A. Plasma physics conditions in E-region ionosphere 

The majority of specular and non-specular radar echoes are 
observed in the lower E region at altitudes 90-120 km, where 
the neutral atmosphere is many orders of magnitudes denser 
than the ionospheric plasma. While the plasma left behind 
a fast meteoroid in the form of a diffusing trail can be sev- 
eral orders of magnitude denser than the background iono- 
spheric plasma, it is usually much less dense than the neu- 
tral atmosphere. The neutral atmosphere remains essentially 
undisturbed by the plasma trail formation. 

All characteristic spatial scales of the diffusion state are 
much larger than the Debye length. The typical diffusion 
time scale is much longer that the electron plasma period. 
Due to these, the quasi-neutrality holds to good accuracy, 
N e w Ni = N, where N e< i are the electron and ion densi- 
ties, respectively. 



In the lower E region, the charged particles collide predom- 
inantly with neutral particles rather than between themselves. 
The following inequalities hold: 57; <C and v en <C O e , 
where Q e> i — eBo/m e ,i are the electron and ion gyrofrequen- 
cies; u en and Vi n are the electron-neutral and ion-neutral colli- 
sion frequencies respectively; e is the elementary charge, m e ,j 
are the electron and ion masses, respectively, and Bq = |Bn| 
where Bo is the geomagnetic flux density. The above inequal- 
ities express the fact that electrons are strongly magnetized, 
while ions are demagnetized due to frequent collisions with 
neutrals. In the equatorial region, B ~ 0.25 x 10~ 4 T, while 
at mid- and high latitudes Bo — 0.5 x 10 _4 T. In the E-region 
ionosphere, 2 + and NO + ions dominate: mj ~ 30 m p , 
where m p is the proton mass, so that rrn/m e ~ 5.5 x 10 4 . 
Throughout the E region, v en ~ 10 Vi n . 

Under conditions Qj <C v in and u en <C fi e > the electron 
and ion diffusion responses to the external electric field differ 
significantly. The response of unmagnetized ions is nearly 
isotropic, while the responses of strongly magnetized elec- 
trons differ dramatically in the directions parallel and perpen- 
dicular to the magnetic field. The ion isotropic mobility gi 
(defined via Vj = <?;F/e, where F is the external force) is 
given by ^ sa e/(mjWj n ). The corresponding electron mo- 
bilities represent tensor components determined via similarly 
defined relations between the electron velocity V e and the ex- 
ternal force F. In the diffusion approximation, the electron 
velocity V e is determined by the inertialess fluid momentum 
equation (in the neutral atmosphere frame of reference), 
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where E is the electric field, P e « NT e is the elec- 
tron pressure, and T e is the electron temperature (in en- 
ergy units) assumed constant. The mobility of magnetized 
electrons along the magnetic field is <jui w e/(m e v en ), 
while in the perpendicular direction the electron mobilities 
e/(m e tt e ) = l/B (the Hall mobility), and 
ev en / (m e Og) (the Pedersen mobility). While the par- 
allel and Pedersen velocities, V e ii and V e j_, have the same di- 
rections as the corresponding components of the electric field, 
the Hall velocity, i.e., the E x Bq drift, is perpendicular to E. 



B. Ambipolar diffusion of plasma columns 

At the earliest stage of trail plasma formation, kinetic pro- 
cesses associated with ionization of ablated material domi- 
nate. After a short time, however, the newly formed plasma 
cools down, typical transport velocities become much smaller 
than the ion-acoustic speed, and the trail spreads over dis- 
tances at least several times the characteristic mean free path 
120]. This can be considered as the beginning of the diffusion 
stage. In this paper, we restrict our treatment to this stage. 

The meteoroid velocity is usually much higher than typical 
diffusion velocities. This means that the trail diffusion starts 
roughly simultaneously over a sufficiently long distance along 
the trail. Considering a part of the trail which is smaller than 
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FIG. 2: Geometry of the trail and magnetic field. 

the typical variation scale of ionospheric and neutral atmo- 
spheric parameters, we will neglect the spatial inhomogene- 
ity and assume approximate translational symmetry along the 
trail axis, making all variations occur in the plane perpendic- 
ular to this axis. 

Assuming a non-zero angle between the trail axis and the 
magnetic field, we will consider the diffusion of a 2D plasma 
trail with the following geometry, see Fig. 13 The trail density 
depends on the x and y coordinates, while it is invariant along 
the z coordinate. The homogeneous magnetic field Bo lies in 
the y, z plane. In this geometry, the plasma density gradient 
V7V and the polarization electric field E = — V$, where <£> is 
the electric potential, have only x and y components, while the 
electron drift velocity may have all three vector components. 

In the special case when the trail axis is strictly parallel to 
B [9], the ambipolar diffusion is axially symmetric around 
the z-axis and its rate, determined by the ratio between the 
ion and electron Pedersen mobilities, roughly follows the low- 
est mobility. The collision frequencies of both electrons and 
ions are proportional to the neutral density which exponen- 
tially decreases with increasing altitude. The ratio of the elec- 
tron Pedersen mobility to the isotropic ion mobility is given 
by the parameter 
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At higher altitudes (usually above 97 km at the equatorial re- 
gion and above 94 km at high latitudes) where the parameter 
ip is less than unity, the lowest is the electron Pedersen mo- 
bility, so that the trail diffusion is determined by the slightly 
increased electron Pedersen diffusion rate. At lower altitudes 
where ip > 1, the lowest is the ion mobility, so that the trail 
diffusion is determined by the slightly increased ion diffusion 
rate. 

In the general case, the ambipolar diffusion is more com- 
plicated because all components of the polarization electric 
field are determined by the same scalar potential, <E>. This 
means that the diffusion in one direction necessarily affects 
the diffusion in other directions. Quasineutrality requires the 
divergences of the electron and ion fluxes to be equal, while 
the fluxes themselves may differ. Nevertheless, vector com- 
ponents of the electron and ion fluxes in any direction usu- 
ally remain comparable. For 2D trail diffusion, if the an- 
gle between the magnetic field and the trail axis is not too 
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small (see the conditions in Sect. Ill C> . there is a signifi- 
cant component of the electron pressure gradient along the 
magnetic field. Accordingly, there should be a significant 
component of the polarization electric field in this direction, 
E\\ = — V||$. The electron mobility along the magnetic field 
Bo (the 'parallel' mobility) is high, while the electron mobili- 
ties in the directions perpendicular to Bo are much lower. The 
electron parallel mobility is much greater than that of ions, 
9e\\l9i ~ miV in /m e v en ~ 5500 (Sect. Ill Ai . Because the 
average fluid velocity of the highly mobile electrons should 
be comparable to that of the low-mobile ions, according to 
Eq. 0, the parallel component of the polarization electric 
field should nearly cancel the electron pressure gradient, 
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where n(x,y) — N/Nq and No is the undisturbed iono- 
spheric background density assumed constant and uniform. 
For isothermal electrons, the electric field and pressure can 
be combined into one force (T e + Ti) V0 rcs , where </> res is a 
dimensionless 'residual' electric potential defined as 



e$ — T e In n 

T P + Ti 



(4) 



Because in our geometry the magnetic field direction has a 
finite y component, while the trail and fields are homoge- 
neous along z, Eqs. l|3} and @ show that the residual potential 
within and around the trail is essentially independent of the y 
coordinate, 



(j> Tes {x,y) « (j> res {x). 



(5) 



The concept of the residual potential cf> rcs is crucial for our 
numerical and analytical treatment. Furthermore, the macro- 
scopic force that drives plasma instabilities is the total force 
acting on electrons, which is just proportional to the gradient 

of (f) Tes . 



C. Restrictions on orientation of trail with B 

The effective cancelation of the electron pressure along the 
magnetic field described by Eq. should only occur if the 
angle between the meteor trail axis and the magnetic field 9 is 
high enough, 
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1.35 x 10~ 2 , (6) 



where the critical angle in degrees is 9o — 0.8°. This means 
that there is a significant component of the electric field along 
Bo to move electrons. The opposite limiting case of 9 <C Oo 
corresponds to the diffusion at nearly perfect alignment be- 
tween the trail axis and magnetic field when no parallel elec- 
tron transport exists 1 10] . However the latter case, as well 
as the most difficult for analysis intermediate case, 9 ~ Go, 
are degenerate cases because of rather small value of the crit- 
ical angle 8q- It is evident that the overwhelming majority 



of meteors penetrating the Earth's atmosphere satisfy the con- 
dition given by Eq. 0. We will also require that the elec- 
tron Pedersen velocity component along the y remains much 
smaller than the corresponding parallel velocity component. 
This yields the condition 
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where the altitude dependent parameter ip was defined in 
Eq. (0. For further simplicity, we will also require 



0> 



(7b) 



The parameter if> decreases exponentially with increasing al- 
titude, see, e.g., Fig. 2 in Ref. ll2~fll . At the altitude range 
where electrons are magnetized while ions are unmagnetized, 
roughly between 80 and 120 km, the parameter if) varies be- 
tween 10~ 3 and 10 4 , so that restrictions given by Eqs. ( fTat and 
(I7bi combined are stronger than (jfjl. In our major analytical 
treatment and simulations, we will discuss the particular case 
of 9 = 90°, but in Appendix 1X1 to find explicit expressions 
for the self-similar solution, we will consider a more general 
case restricted only by Eq. @. 



IIL DIFFUSION EQUATIONS 

If there is no significant ionization and recombination dur- 
ing the diffusion stage, then the line plasma trail density along 
the z axis remains nearly constant. The diffusion stage is ad- 
equately described by fluid equations which include two in- 
ertialess momentum equations for electrons and ions and two 
continuity equations. The quasineutrality condition makes the 
Poisson equation for the electric potential unnecessary. This 
set of equations can be readily reduced to two coupled non- 
linear partial differential equations (PDEs) for the common 
plasma density and the electric potential, 



d t n + V • Ti = 0, 
V • Ti - V • r e = 0, 



(8a) 
(8b) 



where for unmagnetized ions and strongly magnetized elec- 
trons the diffusion flux densities, T e ,i, are given by 



V (e$ + T, t In n) 
n 

V11 (e$-T e km) 



(9a) 
(9b) 



r e _L = n 



i/ ra Vi (e$ - T e In n) bxV± (e$ - T e In n) 



mSl e 



(9c) 



Here the subscripts || and _L pertain to the directions parallel 
and perpendicular to the magnetic field Bo, respectively, and 
b is the unit vector in the Bq direction. The first term in the 
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right-hand side (RHS) of Eq. ( l9cl describes the electron Ped- 
ersen flux, while the second term describes the electron Hall 
flux. Here we assume the isothermal approximation which is 
justified by frequent collisions with the huge thermal reservoir 
of neutral atmosphere (in a more general adiabatic process, we 
would have additional factors 7 6) , in front of T e .% In n.). 

Passing from $ to the residual potential, Eq. l|5}, we rewrite 
Eqs. (|8} in the x, y coordinates as 

d t n - D [V 2 n + V ■ (nV0 res )] = 0, (10a) 

(1 + ip) d x {ndxfces) + Qd v (nd v <j) ICS ) 
+/i (d x <f> ies d y ii - d x n d v cf) ies ) + V 2 n = 0, (10b) 

where we introduced the ambipolar diffusion coefficient, 

T e + T t (T e + Ti) 6 



D- _ . 

a large dimensionless parameter 



(11) 



Q 



sm 



B 2 



sm 2 6» > l,ip, -, (12) 



. ni, i', , ) tj) 
and a dimensionless Hall parameter 

fl= - = — cos6>. (13) 

Here we used the definitions of rp, Eq. (0, and 6o in radians, 
Eq. i|6}. We also used Eq. (fTat and neglected small additions 
to the large parameter Q, which are associated with the ion 
and electron Pedersen mobilities along y. 

For the analytical treatment, it is convenient to rewrite 
Eq. (II Oi in self-similar variables, 

x y 



(Dt) 



1/2 > 



c = 



(Dt) 



1/2- 



(14) 



As a result, we arrive at the following equations for £, t) 

and rcs (£,C,i): 

tdt n - ^ n + ^ n _ v 2 n - V • (nV0 rcs ) = 0, (15a) 

(1 + V) <9 5 (nd£(f) Ies ) + Qd c (nd c <f> rcs ) 

+ [l (<%0 res <9c n - 9g0 res <9 c n) + V 2 n = 0, (15b) 

where the new V operator pertains to the variables defined by 
Eq. d!4i . V = d(). Our solution of meteor trail diffu- 
sion applies this set of coupled PDEs for <^ res and n. Given 
homogeneous background plasma, neutral atmosphere, and 
magnetic field, we assume the following asymptotic bound- 
ary conditions, 



1 



and 







as 



x, y — * ± oo. 



(16) 

We also note that our problem has the following symmetry, 

n(x,y) = n(-x,-y) , (17a) 

0rcs (X, y) = (prcs (~X, -y) . (17b) 

The following section describes a numerical solution of these 
equations and the next section describes their analytic solu- 
tion. 



IV. SIMULATIONS OF TRAIL DIFFUSION AND FIELDS 

In this section, we discuss results of our numerical solu- 
tion of Eq. J 10b using a finite-element PDE solver FlexPDE 
ll3"Tll . The challenge of these simulations is to simultaneously 
resolve both the relatively small scale of the trail density vari- 
ations and the large scale of the residual potential variations 
parallel to the magnetic field (along y). This requires the box 
size along x to be at least several times the effective trail size 
in that direction, a x , while the box size along y should be sev- 
eral times Oq^ 1 , i.e., more than two orders of magnitude, larger 
than the box size along x. FlexPDE uses an adaptive finite ele- 
ments mesh in regions with high gradients to resolve the fields 
and densities with high precision, while uses a coarse mesh in 
regions where gradients remain small. 

We tested the effects of a finite simulation box on the solu- 
tion by varying its size. We also varied boundary conditions, 
setting either the density disturbances, An(t) = n(t) — 1, and 
<pxes, or the corresponding flux densities, Eq. 0, to zero on 
the boundaries. These tests demonstrated that, for sufficiently 
large box sizes, the solution in the inner region remained es- 
sentially unaffected by the choice of boundary conditions. 

In numerical simulations, as well as in the analytical theory 
(Sect. |VJ, we have explored the strictly perpendicular case of 
9 = 90° (p = 0). In this case, we solved Eq. lQj)J with Q 
.").">( )0 corresponding to //;, /»/, = 5.5 x 10 5 and v en /v in = 10, 
see Sect. Ill Al As the initial condition at time t = to, we chose 
a narrow and dense column of plasma, described by the self- 
similar solution (SSS). We used normalized units where the 
diffusion coefficient D = 1 and the initial time for the SSS 
solution to = 1, so that the initial spatial distribution of the 
total normalized plasma density n = N/Nq was given by 



n(x, y, 1) = 1 + Ano exp 



[l + ip)x 2 



2 1 



(18) 



The characteristic sizes of the initial Gaussian density distri- 
bution of the trail in the x and y directions, g x q — [2-0/(1 + 
0)] 1 / 2 and <jyo = V%, are nearly equal for ip I (lower alti- 
tudes), but differ significantly for ip -C 1 (higher altitudes). To 
check the effect of the initial conditions, we tried different ini- 
tial Gaussian density distributions corresponding to the same 
trail line density. We have found that, after the time needed 
for the trail to diffuse over a distance several times the origi- 
nal size of the trail, the solution becomes virtually the same. 
Hence, it is only weakly sensitive to the actual initial peak 
cross-section. 

Equation JlObl involves no time derivatives, so that for- 
mally (j) Ies needs no initial condition. However, the FlexPDE 
application requires setting initial conditions for all variables. 
In our simulations, we usually set up r os(^o) = 0. We tested 
that after a rather short time, the solver automatically sets up 
a time-dependent spatial distribution of <f> Tefj which proves to 
be independent of the initial condition for res . 

Bearing in mind the symmetry along the y direction (i.e., 
parallel to Bo), we simulated a half of the entire space (y < 0) 
with boundary conditions at y = given by zero derivatives 
of both n and <fi lcs . The box size, |y| max = 1000, was al- 
ways at least two orders of magnitude larger than a y (t), but 



FIG. 3: Typical structure of plasma density (a) and residual potential (b). Here ip = 0.2, Ano = 100, and t — 2.5 (the background density 
corresponds to n — 1 and the initial time for the self-similar solution to = 1). 
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FIG. 4: Contours of trail density in the x, y-plane for ij) = 0.05, Ano = 100, t = 5.5: (a): self-similar solution, (b): simulation. 



we needed to keep it that long because of the long-extended 
residual potential. 

In order to model different stages of trail diffusion we per- 
formed several overlapping runs by varying the initial peak 
densities from Ano = 10 4 to Ano = 10- We usually finished 
each run at t = 40, so that the box size x max = 30 remained 
at least several times a x (t). Because of higher anisotropy and 
stronger gradients, the cases of small t/j are much more time 
and resource consuming than those of ij) > 1. On a 3GHz 
Pentium 4 Windows-based PC, our runs lasted from several 
hours (for i/j = 10) to more than a week (for tp = 0.05). Un- 
fortunately, FlexPDE intrinsic restrictions have not allowed us 
to simulate rp < 0.05. 

To monitor the accuracy of simulations, we used general 
relations for the residual potential outlined in Appendix |D] 
and the exact property of density disturbance integrals (DDIs) 
which is analytically derived in AppendixlEl 



Figure [3] shows typical plasma density and residual poten- 
tial plotted over the small fraction of the simulation box near- 
est the trail a short time after the simulation begins. Figure^ 
compares an example of the trail density contours in the x, y- 
plane corresponding to the SSS (a) and to the full simulation 
(b) after the trail has diffused to several of its initial radius. 
Our simulations have revealed the following major features: 

1 . Initially, the plasma density distribution within the trail 
closely follows the anisotropic (for ip < 1) SSS. At 
a later stage, as the peak density falls, the trail re- 
mains nearly Gaussian in each direction but becomes 
more isotropic than does the SSS, see Fig. |4] and it 
diffuses faster, with the isotropic ambipolar diffusion 
rate given described by Eq. il It . The transition from 
anisotropic to nearly isotropic diffusion usually takes 
place while the peak density remains much larger than 
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FIG. 5: Typical variation parallel to the magnetic field (along y) at 
x = of residual potential (solid curve, left scale) and trail density 
(dashed curve, right scale) [for ip — 0.2, Ano = 10 4 , t — 3]. In the 
near-trail region, the density has exponentially strong variation along 
y, while the residual potential varies only within several percent of its 
maximum absolute value. At some point beyond the trail (here about 
ii — — 10V the notential aenrrires a noticeable. 7/-Herivative (electric 




FIG. 6: Disturbances of background density beyond the trail perpen- 
dicular to the magnetic field (along x) [for tjj — 0.2, Ano = 10 4 , 
t — 3.5 at the coordinate \y\ =20 located well beyond the trail]. 



the background density. Section [V] discusses our ana- 
lytical model for this point. 

2. The residual electric potential spreads along y well be- 
yond the trail, see Fig. [5] with spatial gradients along 
y much less than those along x. Along x it has a non- 
monotonic symmetric structure with a deep minimum 
at the trail center and two symmetric maxima around it, 
as shown in Fig. |3lb). At the later stage, when the trail 
diffusion becomes nearly isotropic (see Feature 0, the 
residual electric field becomes much smaller than that 
for the SSS. 

3. Beyond the trail, where the exponentially small trail 
density is much less than the background density, the 
residual electric field extending along Bp may substan- 



tially disturb the density as shown in Fig. |6] The field 
evacuates plasma from the region that maps along Bo 
the potential minimum back to those about the potential 
maxima. As a result, the perturbation density evolves 
in antiphase to the potential distribution along with a 
maximum at x = and two symmetric minima roughly 
at x where the residual potential has the two maxima, 
as visible on the edge of Fig.[3ja). In our simulations, 
the maximum central density compression more than 
doubled the undisturbed background density, while the 
depletions reached about 80% of that. These strong dis- 
turbances were reached near the trail when the peak trail 
density was at least several orders of magnitude larger 
than the background density, see Fig. [5] In all other 
cases or locations, relative density disturbances were 
small. We estimate them in Sect. IV Dl 



We will discuss these features in Sects. IV1 and I VIII when 
describing our analytical theory and comparing it with simu- 
lations. Here we only give brief explanations to some features. 

Feature^ At the later stage of trail diffusion, when the 
initially dense trail density becomes much less dense (but re- 
mains much denser than the background plasma), its gradual 
isotropization and acceleration of diffusion to the ambipolar 
isotropic rate are due to the sharp decrease with time in the 
residual electric field (Feature|2j. In this case, the total force 
acting on electrons becomes so small that diffusion is mainly 
determined by the total pressure gradient acting on ions and, 
hence, is virtually unaffected by the external magnetic field. 
The significant reduction of the residual electric field is due to 
the response of the background plasma which is not included 
in the SSS. To provide quasineutrality, a less dense trail starts 
attracting more charged particles from the background, thus 
reducing the need for strong ambipolar electric field. 

Feature^ The extension of the residual potential along 
y, i.e., parallel to the external magnetic field Bo, is due to 
the high electron mobility along Bo- The deep potential min- 
imum in the trail center is formed by anisotropic ambipolar 
trail diffusion. These two factors are included in the SSS and 
are independent of the existence of background plasma. How- 
ever, the two maxima seen in the residual potential distribution 
along x, see Fig. |3jb), are due to background plasma. The 
electron and ion fluxes originating in the trail extend into the 
background ionosphere with sharply anisotropic and quite dif- 
ferent patterns, while the divergences of the two fluxes remain 
balanced. Beyond the trail, these patterns have a quadrupole- 
like structures providing current closure. To drive the return 
currents, the background ionosphere develops potential gradi- 
ents which oppose those within the trail, i.e., those responsible 
for the trail diffusion. This gives rise to the two symmetric po- 
tential bumps around the deep potential minimum and draws 
background plasma into the trail edges. 

The simulations allow us to determine diffusion rates and 
isotropy for a range of altitude-dependent ijj. This will be pre- 
sented in conjunction with our analytical theory in Sect. lVIll 
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V. ANALYTICAL THEORY 

Here we discuss our analytical theory for the 'perpendicu- 
lar' case of 8 — 90° (/i = 0), when a double mirror symmetry, 



n(£,C)=n(±£,±0, 



(19a) 
(19b) 



takes place. The Hall velocity of electrons is directed strictly 
along the trail axis, so that it does not affect the 2D diffusion 
of the plasma trail. In the general case of 8 ^ 90°, the Hall 
mobility breaks the double mirror symmetry but keeps a rota- 
tional (by 180°) symmetry around the trail axis, Eq. dl7t . For 
sufficiently large angles 8, distortions caused by the Hall mo- 
bility are relatively small (see AppendixlATi. so that the double 
mirror symmetry roughly takes place. However, the Hall mo- 
bility affects the aspect ratio of typical scales along each axis. 
We reserve the more general case for future work. 

The trail diffusion and evolution of the ambipolar fields are 
described by two coupled nonlinear PDEs for the plasma den- 
sity and the residual potential, see Eqs. i ll Ob or d!5l >. Our the- 
ory provides an approximate analytical solution of both inter- 
related problems. The key point of the theory has been the 
inclusion of the background plasma into consideration. In so 
doing, even for p, = 0, one cannot find the exact analytical 
solution of the coupled equations. However, based on the ex- 
istence of the large parameter Q — 8q~ 2 ~ 5500, Eq. dl 2L 
and on the insight from our simulations, see Sect. llVI we have 
developed an approximate approach based on the perturbation 
technique. 

We have used the fact that the initially narrow and dense 
trail keeps nearly Gaussian density distribution during the trail 
lifetime even if the diffusion rate is changing with time, see 
Sect. IIVI Given the characteristic scale of the Gaussian dis- 
tribution, we have reduced Eq. < l 1 5bl to a simpler equation 
for the residual potential. This one-dimensional (ID) linear 
integrodifferential equation is the governing equation of our 
theory. It has two equivalent forms that depend upon only one 
dimensionless parameter p. This master parameter is propor- 
tional to the the square of the characteristic Gaussian scale 
along x, p ix. cr^{t), see Eq. ( I65t below, and, hence, monoton- 
ically grows with time. Approximate solutions of the govern- 
ing equation in different domains of p(t) have allowed us to 
obtain the spatial distribution of the residual electric potential 
in the entire 2D space. 

The trail diffusion is mainly determined by the parabolic 
behavior of the potential along x within the central region and 
is practically insensitive to the potential behavior outside this 
region. This has allowed us to find from Eq. dl 5a> the ex- 
plicit time dependence of p(t) and thus close the solution. As 
a result, we have obtained approximate analytic expressions 
describing the trail diffusion and the evolution of the residual 
electric fields. These expressions reasonably well agree with 
the numerical results and can be used for quantitative predic- 
tions. 

Our analytical theory consists of several major steps: (1) 
obtaining explicit expressions for the self-similar solution; (2) 



solving an equation for the electric potential in the 'far zone' 
where the density disturbances are relatively small; (3) ob- 
taining the governing equation for the residual potential in the 
'near zone' depending on the master parameter p; (4) solution 
of this equation in different ranges of p; (5) finding expres- 
sions for p(t) in the initial value problem, which gives the 
approximate solution of the trail diffusion problem and closes 
the solution for the residual potential; (6) estimating density 
disturbances of the background plasma beyond the trail and 
obtaining a correction factor for p ^> 1 at a later stage of 
trail diffusion. We have implemented steps (1) to (3) in Ap- 
pendices |X| to |Cj while in Sect. IV Al we will only explain the 
basic ideas and discuss the resultant governing equation. We 
have implemented steps (4) to (6) in Sect. IV Bl to lV Dl 



A. Governing equation for near-zone residual potential 

Here we briefly describe the governing equation intended 
for the solution of the residual potential problem. Given the 
trail density spatial distribution, this approximate equation is 
derived from Eq. dl5bl i. The derivation requires separation of 
the entire £, £ plane into two overlapping regions: the near 
zone, 



ICI 



\v\ 



< Q 



1/2 



and the far zone, 



ICI » i. 



(20) 



(21) 



We emphasize that the terms 'near' and 'far' pertain to the 
coordinate £ only, so that the near zone includes the entire £- 
axis (the rc-axis). The residual potential in the near zone is £ 
dependent but remains nearly constant in the £ direction, i.e., 
along B , see Eq. (|5}, making its calculation a ID problem. 

In the far zone, which includes the background ionosphere 
only, the residual potential has a slow £ dependence, but we 
may roughly neglect there density disturbances (see Appendix 
[B]and Sect. I V DV The fact that divergences of the electron and 
ion fluxes are nearly equal, as required by quasi-neutrality, 
results in the Laplace equation in terms of some renormal- 
ized coordinates. In those renormalized coordinates, the en- 
tire near zone reduces to a thin cut which provides boundary 
conditions for the Laplace equation in each of the two half- 
spaces, £ > or £ < 0. Because of the symmetry, it is suffi- 
cient to consider only one of these half-spaces (£ > 0). Given 
the residual potential in the near zone 0j? cs , a straightforward 
solution of the Laplace equation in integral form yields the 
spatial distribution of the residual potential in the entire 2D 
space, as described by Eq. iB4\ . 

Using the Gaussian approximation of the trail plasma den- 
sity, Eq. dC7t . 
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(22) 



and integrating it over the near zone, we obtain the sought-for 
boundary condition for the Laplace equation. This integration 
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involves 0° cs (with still unknown £ dependence) multiplied 
by the corresponding trail density integral over the coordinate 
£. Associating the upper half-space with a complex plane in 
terms of renormalized coordinates and using analytical prop- 
erties of the potential in the far zone (see Appendix iBi. in 
Appendixlclwe obtain for the the near-zone residual potential 
an integrodifferential equation in two different, but equivalent, 
forms: 



2p_ 

71- 3 / 2 



<p(r) 



dr + d v (p{rf) = 77, 



(23a) 



p<p(v) - 



2y^ 



dMT) .e-^dr = S(r,). (23b) 



r — ?] 



Here the renormalized near-zone residual potential <p and co- 
ordinate r\ are defined by 



pfa) = 2(1 + ^.(0, 



(24a) 



V 



£ x 

/V 5 = (qDt) 1 / 2 



'-!=■. (24b) 



The coupled dimensionless parameters q(t) and p(t) are pro- 
portional to the square of the Gaussian peak dispersion with 
respect to x coordinate, a 2 , Eq. d51bl . 



p(t) = iq {t)t 



IK 

2D 



where 
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(25) 



(26) 



and Ny m is the line plasma density of the trail along its axis z. 

The important function S(rf) = rje^ r ' ^C!^ e T dr — 
1 in the RHS of Eq. d23b1 . see also Eq. JC171 together 
with the highest-order power-series approximations of S(rj), 
Eqs. JC20l . OH . is shown in Fig. El In Sect. lVBl we will 
see that this function is a good qualitative, and in some cases 
even quantitative, representation of the residual potential in 
the near zone. 

Equivalent Eqs. ( I23b .b) are linear integrodifferential equa- 
tions in renormalized variables defined by Eq. (1241 . In 
these variables, these equations depend upon only one di- 
mensionless parameter p. We will use one of two equivalent 
Eqs. (I23h .b) depending upon the convenience for the specific 
analysis. In Section|D] we discuss the general analytical prop- 
erties of the solution, which can be used for monitoring the 
solution. In particular, under assumed approximations, the in- 
tegral of <^>rcs(£) over £ should be equal to zero. This explains 
mathematically the non-mono tonic, two-bump structure of the 
residual potential in the near zone observed in our numerical 
calculations (see Sect. II VI Feature 2). 

The physical sense of the outlined mathematical procedure 
is that different patterns of the electron and ion fluxes in the far 
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FIG. 7: Function S{>q) = r/e" " 2/4 £ /2 e^dr - 1, Eq. JcTTi (solid 
curve), along with its parabolic approximation, — 1 +77 /2, Eq. <C20t 
(dash-dotted curve), and large-77 asymptotics, 2 /if, Eq. <C2H 
(dashed curve). 



zone create a feedback to the near zone. The non-monotonic 
structure of the residual potential (in each of the two half- 
spaces, x > and x < 0), can be explained as follows. The 
total force acting on negatively charged electrons is propor- 
tional to eV0° os . Within the major trail, this force via elec- 
tron Pedersen diffusion pushes electrons outwards. This re- 
quires the residual potential to have a minimum in the trail 
center. Well beyond the trail, different electron and ion flux 
patterns form return fluxes in the background ionosphere. On 
approaching the trail, these fluxes are directed inwards which 
requires oppositely directed gradients of the residual potential. 
This gives rise to the two bumps of the residual potential. 



B. Solution of governing equation 

Having solved equations for the residual potential in the 
near zone and using formulas of AppendixlHl we can find the 
ambipolar electric field and estimate the plasma density dis- 
turbances both in the trail and in the background ionosphere. 
Equivalent Eqs. d23fa.b) have singular Cauchy kernels. While 
gen eral methods for solving singular integral equations exist 
l22l I23I1 . we are unaware of such methods for Cauchy-type 
integrodifferential equations. Even the existence and unique- 
ness of the solutions of such equations is not a priori known 
and should be studied individually for each specific equation. 
Nevertheless, for each value of p, Eq. d23l > has the unique so- 
lution. We will not dwell here on the proof, but note that the 
uniqueness of the solution is provided by the positive sign of 
p [for negative p, Eq. d23i would have infinite number of so- 
lutions]. 

While the equation for 92(77) has the unique solution, it is 
not possible to find it analytically in the general case. In this 
section, we discuss approximate solutions of Eq. d23l l for var- 
ious ranges of p. We will proceed from the simplest case to 
the more complicated ones. 
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1. Late stage diffusion, p> 1 

The simplest limiting case of p ^> 1 corresponds to a 
sufficiently late stage of the meteor trail diffusion, when the 
peak trail density exceeds the background not too signifi- 
cantly, (iV max - N )/N « ^/ 2 /[2(l + ^)6 ] (V 1/2 /©o = 
V in/Qi ^> 1). In this limiting case, one can solve Eq. (123b) 
using a perturbation technique. For p 1, the integral term 
in Eq. J23bi is small compared to the first term in the RHS. To 
the zeroth order, we can neglect the integral term so that the 
zero-order solution is given by ip(rj) (p(°> (rj), 




^ 0) (V) = (27) 
P 

In the central region, \r}\ <C 1, the renormalized amplitude of 
the electric field, d v ip ~ rj 2 /(2p), proves to be much less than 
that for the self-similar solution, d v ip ss = rj 2 /2. 

The physical reason for much smaller residual electric field 
is as follows. When electrons leave a sufficiently dense trail, 
mainly along the magnetic field, slowly diffusing ions create 
a significant ambipolar electric field to retard electrons. If the 
trail is not dense then there are enough background electrons 
to substitute for those leaving the trail, so that no strong resid- 
ual electric field is needed. The integral term in Eq. J23bb . 
corresponding to the term d v ip in Eq. d23al >. describes the elec- 
tric field associated with the trail electrons, while the remain- 
ing terms in the left-hand sides (LHSs) include the effect of 
background electrons. In the case of p ^> 1, the latter terms 
dominate. 

To verify that the neglected integral term is really small, we 
need to obtain the next-order approximation. To this end, we 
substitute the zeroth-order expression tp^ (77) to the integral 
term in Eq. J23bi . Using Eqs. dC22t -( IC24l . we obtain the 
first-order approximation, 

^)(,) = «- (28) 
P P P P 

where J(rj) is given by Eq. dC23> and I ps 0.643, p ~ 1.546, 
Eq. dC25> . Figure [8] shows the function J(rj) and its approx- 
imate counterpart J(rj) = lS(prj). The two functions prac- 
tically perfectly agree everywhere except the two symmetric 
maxima. While the approximate expression in the RHS of 
Eq. d28l is less accurate, it is much simpler for analysis. Equa- 
tion d28l clearly shows that since the second term is small 
compared to the first one the perturbation approach employed 
here is consistent. 

As seen from Eq. d28i . for p ^> 1 the residual potential is 
approximately a linear combination of two S'-functions with 
different arguments. The function S(rj) has major qualita- 
tive properties of the residual potential shown in Fig. 0t>). 
One can also check that this function automatically satis- 
fies the general analytical properties of the solution given by 
Eqs. l|DT}, «|D4j, and (fOTl. 



FIG. 8: First-order corrections in tp^ (rj), Eqs. 1281 and 13 1 i . mul- 
tiplied by — p 2 . Solid curve: J(rj), Eq. IC23I : dashed curve: J = 
lS(prj), Eq. <C24l : dot-dashed curve: the second term (multiplied by 
-p 2 ) in the RHS of Eq. f or o » 1. 



2. Intermediate case, p ~ 1 

The perturbation technique of the previous subsection fails 
for p ~ 1, and especially for p -C 1. In these cases, we have 
found no rigorous analytical solution. However, we have de- 
veloped an approximate approach for solving the integrodif- 
ferential equations using analytical fitting with iterations. This 
approach is similar to the numerical shooting method. Assum- 
ing a reasonable initial approximation for a zero-order itera- 
tion ip(rj) ~ <^(°) (rj) with adjustable parameters, we may put 
this function into the integral term of one of the integrodiffer- 
ential equations and find the next-order iteration. Comparison 
between the two iterations at the critical central region allows 
us to adjust the parameters in order to find the best fit. 

The two forms of the integrodifferential equation for the 
near-zone potential, Eqs. (I23h .b). are equivalent in rigorous 
mathematical sense, but this is not so if we seek the trial- 
and-error approximate solutions. For the intermediate case 
of p ~ 1, it is more convenient to use Eq. d23bl . This allows 
us to disregard the correct asymptotic behavior, Eq. SD41 . of 
the initial trial function [this behavior is required in Eq. d23a> 
to eliminate the diverging effect of the exponentially growing 
factor in front of the integral] . 

We start by choosing a simple trial function, ip^ (77), which 
satisfies the condition (ID1> but not necessarily dD4t. Putting 
it into the integrodifferential term in Eq. d23bi . we calculate 
the first-order iteration, 

p V 2 V7r y-r) J 

(29) 

It can be readily shown that ^^^(77) will automatically sat- 
isfy Eqs. dD li and JD4> . In principle, the same procedure 
could be continued further to calculate next-order approxima- 
tions ip^(ri), (p( 3 \r)), but because of increasing complexity 
the analytical calculation of the integrals is hardly possible. 

To make calculations simpler, we choose as the initial 
function ip(°\r]) the even 77-derivative of the odd function 



11 



-A?7exp(— A77 2 ), 



^ (0) (77) = A (2\rf - 1) exp (-At/ 2 ) 



(30) 



Because the generating function — At] exp (— A77 2 ) tends to 
zero as 77 — > 00, the trial function with <j)(rj) = V'^H 7 ?) au ~ 
tomatically satisfies Eq. (ID1I . Besides, in the central part it 
qualitatively looks like the expected residual potential around 
the major trail. 

Now we find ip^ (77), Eq. (I29> . Direct calculation yields 



4A 
1 + 4A 



2XA 



py/1 + 4A 

+ (3 - 2Ar/ 2 ) S , ( ? 7\/rT4A) 



(31) 



The major residual potential is located in the central region, 
where it is nearly parabolic. It is natural then to fit the 
parabolic behavior of the two functions, (p(°\0) and (p^\ at 
small |ry|. For p > 0.1, this will provide a good fit in the en- 
tire central region and, due to Eq. (ID 111 , a reasonable fit in the 
nearby region as well. Fitting the two parabolas, 

^°)(0) = ^)(0), 9;/)(0) = 3,> (1) (0), (32) 

we obtain two equations for the two unknown parameters, A 
and A. Using the Taylor expansion for the function S(rj) at 
small 77, Eq. dC20t . and similar for y>(°) (77), Eq. POL we ob- 
tain 



A = 



8A 2 (5 + 16A) 
(1-6A)(1 + 4A) 3/2 ' 

(1 - 6A) (1 + 4A) 3/2 
2A(3 + 10A+ 16A 2 )' 



(33a) 
(33b) 



These expressions give an implicit dependence of the two fit- 
ting parameters, A and A, on p. The parameter A varies in the 
range from (p — > 0, A — > 00) to 1/6 (p — ► 00, A — ► 0). We 
cannot explicitly express parameters A and A in terms of p. 
However, we can determine the asymptotic behavior of A(p) 
and A(p) at large and small values of p and construct interpo- 
lating formulas that would provide a reasonable approxima- 
tion to the exact solution. We choose the following interpola- 
tion formulas, 



A ~ A in t(p) 



A« AntO) = 



p + cp- 



1/2 



40 + bp + dp 2 

10 + Np 
9p + Gp 2 + Ap 3 



1/2 



(34a) 
(34b) 



where b Rj 122.7, c « 1.52, d Rj 54.74, A Rj 14.67, and G « 
29.7. These approximations deviate from the exact functions 
A(p) and A(p) only several percent at worst (when p is small). 

For large p, the first-order solution (77) matches the cor- 
responding solution obtained in the previous subsection. In- 
deed, the major terms in Eqs. J2 8t and ( 13 1> are the same, 
S{rj)/p. Furthermore, for p > 1 we have Awl/6 and 



Aw X/p, In this case, the smaller second term in the RHS 
of Eq. d3 1 1 approximately agrees with the second term in 
Eq. d28l with the largest mismatch near the two maxima of 
ip^(r]), \rj\ w 2. Figure |8] shows that the relative mismatch 



there is about 10% for exact ip^ (77) in Eq. d28l > and is roughly 
twice as large for the approximate expression [in terms of 
S(pi])]. In all other locations, the agreement is much bet- 
ter. Furthermore, because this is the mismatch between the 
two minor terms, the relative mismatch between the two full 
solutions ip^ (77), Eqs. d28l and d3 1 1 , is p times smaller. 

Figure |9] shows the fitting solution for p = 1 and p = 0.1 
and the parabolic SSS for the residual potential, <p ss {rj) = 
ip(0) + tj 2 /2. From Fig.|9ja) we see that for p = 1 the two 
approximations, cp^(rf) and <p^'(rf), are reasonably close to 
each other in the entire central region, but the self-similar so- 
lution deviates significantly from both. Figure|9jb) shows that 
for small p the two approximations start deviating from each 
other even in the central region. The reason for this is that 
the lowest parabolic term oc ?7 2 in the Taylor expansion for 
ip^irf) becomes so small that the higher-order terms cx ?7 4 
start playing the major role, even at sufficiently small 77. At 
the same time, the first-order approximation ip^'(rf) closely 
approaches the self-similar solution in the central region, as is 
should for small p (see the following Section). 



3. Early stage diffusion, p <C 1 

The limiting case of p <C 1 deserves special attention be- 
cause it usually applies when the meteor trail is first detected. 
This case corresponds to the early diffusion stage of a suf- 
ficiently dense plasma trail, when the trail density evolution 
follows the SSS, so that q w q = ip/(l + tp), see Eq. ( Bob 
below. According to Eqs. d25i and d26i . if p -C 1 then 



An 



1/2 



A 2(1 + V)6 2(l+i>)Sk' 

This case is the most difficult case for the analysis. Being un- 
able to obtain a rigorous solution in the entire range of ip {rj), 
we can construct a reasonable approximation. Qualitatively, 
the solution has the same basic form as for p > 1, viz., a 
deep potential minimum surrounded by two symmetric max- 
ima. Smallness of p suggests neglecting the integral term in 
the RHS of Eq. j23al , so that we obtain the approximate ex- 
pression 



p (77) Rj p (0) + ^ 



(35) 



corresponding to the parabolic SSS with the unknown <p(0). 
However, this is only valid within a restricted range of 77, 
where the exponentially growing factor in front of the integral 
in the LHS of Eq. (|55aJ is of no importance. As pexp(?7 2 /4) 
becomes of order unity, i.e., the value of |?7| approaches a crit- 
ical value 770 — 2 yTxt (1/p), the fast growing exponential fac- 
tor starts overpowering the small parameter p, so that for suf- 
ficiently large |?7|, the integral term in the LHS of Eq. ( 12 3 al l 
becomes more important. The critical coordinates \rj\ = 7/0 



FIG. 9: Zero-order, <p (t)) {rj), and first-order, p {1) (r)), approximations: (a) p = 1 (A = 0.107, A = 0.675); (b) p = 0.1 (A = 4.52 x 10~ 2 , 
A = 2.97). The dashed parabolic curves in the middle show the self-similar solution, <p BB (r]) — tp(0) + rf /2. 



are located slightly beyond the two maxima of the residual po- 
tential, where the parabolically growing potential of the major 
trail starts transforming to a decreasing potential beyond the 
trail. 

The parameter 770 is of major interest for us because it de- 
termines the rough boundary for the parabolic potential and 
hence for the maximum residual electric field responsible for 
the electron drift. Because the above heuristic estimate for 
7/0 was based on a simple order-of-magnitude comparison, we 
need a more accurate and consistent estimate. We will con- 
struct a piece-wise approximation for <f{rf), which would be 
close to the self-similar approximation given by Eq. d35l at 
|ry| < ?/o an d proportional to 7/~ 2 , Eq. dC2 ll >. at \r)\ > r/Q. 
We will do this by iterations. As a zero-order approximation, 
ip( \r]), we will set 



where 



i( v ) 



(2 2 

[Vo ~V 
2% {Vo 



V 2 + 2<^o) 



In 



Vo + V 



Vo~V 



(39) 



The function 6(pW(rj) is negative with the roughly exponen- 
tially growing absolute value. At \rj\ = 770, we have 



drj 



n=vo 



a 2sl 
tt3/2 ' 



(40) 



which allows us to match the derivative of tpW at both sides 
of M = Vo, 



v {0 Hv) = 




£ ^ if \n\ < vo, 

J f if \ V \ > ryo, 



(36) 



with an unknown constant tpQ. The function <p^(rf) is a con- 
tinuous function, while its derivative, d^tp^ (rj), is discontin- 
uous at 1 77 1 — r/Q. As a next-order approximation, we will 
construct a smooth function tp^- 1 ' (77), 



ipo+^+S^Hv) 2 if \v\<Vo, 
(^o + #+<M 1} M) $ if \v\ > Vo, 



(37) 

where the deviation from the self-similar solution, dip^ (r/), 
can be found from Eq. JD6a> equivalent to Eq. d23al . Substi- 
tuting <^ ' (t) for tp(r) in the integral term, we obtain 



drj 



2pe- 



r 3/2 



dr 



pe"* 1(77) 



7r 3 / 2 



(38) 



Vo + 



4pe 4 tp 



2 
Vo 



<Po 



9 
Vo 



7T 3 / 2 

Integrating Eq. J38l> . we obtain 

P 



5<pM(ri) = 



-3/2 



<V (1) (vo) 



'I(r)dr, 



(41) 



(42) 



where, notwithstanding the singular denominators in I(r]), the 
integrand at r — > behaves regularly, I(t) oc t. Assuming 
sufficiently large 775% to the first-order accuracy in the t]q ex- 
pansion, we obtain 



<V (1) (%) 



r 3/2 



Vo 



1 



6 — 27 — 4 In 770 



Vo 



'•po 



+ 4-27-4111770], 



(43) 



where 7 w 0.5772 is the Euler constant. Substituting this 
to Eq. i41i . we obtain the relation between p, 770, and (po. 
To exclude the unknown constant ipo and determine p as a 
function of 7/ , we need one more equation. We will use the 
general relation ip (770) drj = 0, Eq. (IDli . which is neces- 
sary for consistency of the Laplace equation in the far zone, 
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see Appendix 151 Integrating the major parabolic part of the 
potential at \rj\ < rj in Eq. d37l is trivial, but integrating 
Stpt 1 ) (77), Eq. (1421 . is complicated. Note, however, that in- 
tegrating the major part gives rise to an additional large factor 
~ 770 an d every next integration of Sip^ 1 ' (77) gives rise to an 
additional reduction factor ~ 77^" 2 . While the derivative of 
Sip^ at |?7| = 770, Eq. J38l >. is comparable to 770, the relative 
contributions of Sip^ (r]o), and especially of its integral, are 
small compared to those from the major part, bearing this 
in mind, we obtain to the first-order accuracy with respect to 

Vo 2 , 



p « 



(44a) 
(44b) 



Application of other general relations, Eq. ( ID7L instead of 
Eq. JD li . yields somewhat different values of the first-order 
corrections. In Eq. ( 144 at . we will disregard these corrections 
and obtain the solution of 770 (p) in terms of the Lambert W- 
function, W(x) 01: 



2W 



2p 2 



1/2 



(45) 



To logarithmic accuracy, we obtain 



Vo - 2 < In 



r 3/2 



2pJln(nVyp) 



1/2 



2 In • 



1/2 



in good agreement with the above heuristic estimate. 



the density disturbance An(rj ) 
simulations, 



n(to) — 1, as we did in our 



An(t ) 



n ss (t ) 
An 



Auq exp 



1 fl 

4 



47rD/j V 4> 



i + 4' 



4, 
1/2 



(48) 



Our numerical computations, starting from the initial con- 
dition 7T,Traii (io) = nSS (to) with Ho 3> 1, show that, for 
some time, the trail density approximately follows the SSS. 
As the peak density decreases with increasing time, the solu- 
tion starts to gradually deviate from the SSS. However, within 
the trail density peak, it keeps the nearly Gaussian form, 



"-Trail 



An Q t Q 

t 



go 

?(*) 



exp 



(49) 



where the characteristic diffusion scale in the x direction (£ = 
x/^/lji) is determined by a time-dependent coefficient q(t) 
with the SSS initial condition, 



q(t ) 



(50) 



In the original coordinates x and y, the nearly Gaussian peak 
can be expressed in terms of the x, y-dispersions, <T x ,y< 



cr x {t) 



(x 2 (t)) 

x 2 n Trili i(x,0,t)dx 



(46) a y (t) 



(y 2 (t)) 

y 2 nTva.n(0,y,t)dy 



(51a) 

^Traii(^,0,t)da;, 
(51b) 

nTrai\(0,y,t)dy, 



C. Trail diffusion 

Given the approximate expressions for the residual poten- 
tial, we will solve now the problem of trail diffusion described 
by Eqs. dl0al > or d 1 5 al > . While this treatment will require a 
number of approximations, the comparison with simulations 
described in Sect. IVIll shows that the analytical theory devel- 
oped below is valid to a good accuracy. 

For the strictly perpendicular case, 9 = 90°, the self-similar 
solution (SSS) obtained in AppendixlAlis given by 



as 



7V SS 




= 








x exp 


[-K 1 



1/2 



4> 



i + 4> 

4> 



(47) 



where nii n = N\[ n /No is the trail line density along the axis z 
normalized to the background plasma density. Under assump- 
tion that the initial maximum density is well above the back- 
ground plasma density, N max ^> N , the self-similar Gaus- 
sian profile can be used as a reasonable initial condition for 



"Trail = 



"lir 



2-KU x (J y 



exp 



X 

2V! 



(52) 



The small residual electric field in the y direction (i.e., along 
Bo) practically does not affect the SSS, cr 2 = 2Dt, because 
the parallel residual electric field is small, so that the diffu- 
sion along Bo is determined by the isotropic diffusion rate 
D, Eq. il It . At the same time, the strong residual electric 
field along x affects the diffusion in this direction, especially 
at altitudes above 93-97 (ip < 1). The deviation of the char- 
acteristic trail size from that determined by isotropic diffusion 
is described in Eq. (I49i by the coefficient 



<?(*) = 



2Dt ' 



(53) 



The coefficient q(t) slowly grows with time as the trail dif- 
fuses faster than it does in the SSS. As the electric field falls 
well below that predicted by the SSS (see Sec. lVBt . the dif- 
fusion becomes nearly isotropic, q(t) ~ 1, determined in the 
two directions roughly by the same isotropic diffusion rate D. 
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The coefficient q(t) is determined by the residual electric 
potential in the near zone. In the central part of the near zone, 
the potential along £, </>res(£)> always has a nearly parabolic 
behavior similar to the SSS, Eq. (IA2> . 



res 



c 



(o) 



(54) 



Unlike Eq. iA6\ . however, the coefficient B xx in Eq. (I54i is 
no longer a constant but changes with time. Our simulations 
show that the central region described by the nearly parabolic 
^-dependence is broad enough to include the entire trail peak. 
It is the parabolic behavior of the potential that imposes the 
Gaussian shape of the trail. 

Substituting Eqs. J49> and \5A\ into the continuity equation 
dl5al >. we obtain a differential equation for q(t), 



dq 
~dt 



(1+B X 



(55) 



which is not yet closed because the parabolic coefficient B xx 
is itself an unknown function of q{t). To determine it, we will 
use the solutions for the residual potential found in Sect. IV Bl 
These solutions are expressed in terms of renormalized vari- 
ables <p and r\ defined by Eqs. J24i . In these variables, the 
parabolic behavior of Eq. J54i is described by 



<p(v) ~ <f(ty 



Frf 



(56) 



Here the coefficient F(p) is similar to B xx . Unlike the latter, it 
depends upon one parameter p related to q by Eq. J25l >. pit) = 
jq(t)t. Comparing Eqs. i54l and d56i using i24l yields the 
relation B xx q = F(p)/(1 + tp), which reduces Eq. (1551 to a 
closed nonlinear differential equation for pit), 



1 dp 
7 ~dl 



F( P ) 



+ 1. 



(57) 



Integrating Eq. ( I57t for t (p) with the initial condition at t = to 
given by p(t ) w p = jq t , 



Po 



7^o 



1 + ip 2(l + i/>)9 n 
Eq. ( IC141 . we obtain for the general F(p): 

dg 



< 1, 



7(<-<o) 



i-FG?)/(i + i/0 



(58) 



(59) 



During the initial time when p(t) <C 1, both the trail density 
and the parabolic approximation of the potential follow the 
SSS, B xx » gwft.Fwl, 



P(t) 



jipt 
1 + ip' 



(60) 



Using this, in Eq. ( 1591 we can set F = 1 for t < to and then 
replace to and p by zeroes. This corresponds to diffusion 
of the initially infinitely dense and thin ((5-function like) trail. 
The solution of p(t) is determined by reversing Eq. ( 1591 . The 



RHS of Eq. j59t is always positive, so that this reversal is 
unambiguous. 

In the course of the trail diffusion, when the parameter p be- 
comes of order unity and greater, the monotonically decreas- 
ing function F(p) starts deviating from unity. To obtain the 
explicit expression for F(p) = {l/rj){d<p/drj), we will use 
the solutions found in see Sect. IV Bl While we have no sim- 
ple analytical expression for F(p) for all p, we can approxi- 
mate it by interpolating between F m 1 for p <C 1, Eq. ( 1351 . 
and F w l/p at p ^> 1, Eq. i27\ . The simplest monotonic 
interpolation is 



Ho) 



i 



i + f? 



(61) 



Substituting this into Eq. ( I59l l allows us to obtain the solution 

ioit(p): 



1 , / 1 + tp 
jt = p+ In I H ■ — p 



1 + 



1> 



(62) 



Reversing this relation, we obtain an explicit expression for 
pit) in terms of the Lambert VK-function, 



Pit) 



(63) 



This general expression describes the transition from the SSS 
for jt < 1, Eq. to 



pit) « 7* - 



1 .lnfl±^ 



i + -0 V i> 

at large jt. According to Eqs. ( I25> . ( I51bl i, and ( I53> . 



(64) 



= — I — = Z77 PW> (65) 



7 



7TiV 



so that the first term on the RHSs of Eqs. ( I64i corresponds to 
diffusion over x with the isotropic rate D, while the second 
logarithmic term describes a small time delay. Since the trail 
diffusion along y has the same rate, w 2Di, the above 
expressions predict an evolution from anisotropic to isotropic 
diffusion. 

The diffusion starts changing its character from the self- 
similar, sharply anisotropic diffusion to the nearly isotropic 
one roughly at a time t = t CT when p(< cr ) = 1. According to 
Eqs. d26t and ( I62K the corresponding critical time t CI is given 
by 



trx — 



where 



N Xin eB Q K,{^) 



1000K 



2tt (T e + Ti) 



10 14 irr 



No 



T e + Ti 



Bo \ 

10 4 nT ) 

fC(ip), (66) 




1 + + In 



1 + 2V> S 



(67) 
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The critical transition time from anisotropic to isotropic diffu- 
sion is mostly relevant for high altitudes with ip <C 1, where 
the anisotropy of diffusion is clearly pronounced. In this case, 
the function tCOip) « + ln(l/^)]. According to Eq. (j2ji, 
y/lj) oc NnB^ 1 . From Eq. d66t we see that i cr is practically 
independent of the magnetic field, making the high and low 
latitude diffusion evolve similarly. The altitudinal dependence 
of t CI , however, depends strongly on the background plasma 
density No and even more so on the neutral density N n . 



D. Density disturbances beyond the trail 

The residual electric potential originated within the trail and 
extended along the magnetic density creates disturbances in 
the background plasma, as we observed in our numerical cal- 
culations, see Sect. II VI Feature 3. These density disturbances 
occur due to the fact that the extended residual potential at- 
tracts ions from the surrounding ionosphere into the central 
part with the deep potential minimum around x — 0. This 
results in plasma compression in the central potential mini- 
mum and depletion in the adjacent regions where the residual 
potential has two symmetric maxima. 

In this section, we will estimate density disturbances in the 
background plasma beyond the trail. One motivation for this 
is that when obtaining the equation for the residual potential 
in the near zone, Eq. ( I23K we have completely neglected these 
disturbances. This has allowed us to reduce the equation for 
the residual potential in the far region to the Laplace equation 
in renormalized coordinates, see Eq. (IB 1 1 in Appendix 151 To 
estimate the effect of density disturbances, we will consider 
the strongest case of dense trail described in Sect. IV B 31 An- 
other motivation is to find a proper adjustment for the parame- 
ter p, which is needed for the residual potential at sufficiently 
large state of trail diffusion. 



1. Strongest case, p <C 1.' justification of analytical approach 

To justify the neglect of density disturbances beyond the 
trail, see Appendix 151 we will make a simple analytical esti- 
mate of density disturbances beyond the trail. The largest dis- 
turbances beyond the trail occur in the near zone in the early 
stage of dense-trail diffusion, p <^1. 

Given the spatial distribution of the residual potential 
<^ res , the dynamics of density perturbations is described by 
Eq. dlOal If we assumed a stationary regime, dpn = 0, then 
imposing boundary conditions at infinity where n — > 1 and 
(j)® es — > we would obtain the local Boltzmann distribution 
of plasma density, n — exp(— 0° cs ). In the strongest case 
of p <C 1, according to Eq. d44b> . the potential minimum is 
ipo ~ — 7?q/3, so that according to Eq. ( I24a> the correspond- 
ing potential <^ es « 0° cs = tpo/2(l + ip) ~ -»)g/[6(l + V)]- 
For sufficiently small p, the parameter 770 depends logarithmi- 
cally on p and may reach several units, Eq. ( 146 \ . In this case, 
the local Boltzmann distribution would yield disturbances of 
the background plasma near the trail much greater than the 



background density itself, n ~ 1/p 1, making the as- 
sumption underlying our theoretical approach to be invalid. 
Fortunately, our numerical simulations, even in the cases as 
strong as n max = 10 4 , have shown that the background den- 
sity disturbances are at worst of order unity in the central near 
zone and are much less beyond it, as described in Sect. I1VI 
The physical reason for this is that the trail diffusion is a non- 
stationary process with a diffusion prehistory. Plasma den- 
sity disturbances caused by the potential expanding with time 
prove to be noticeably smaller than those caused by the sta- 
tionary potential. 

For a simple estimate of plasma density disturbances out- 
side the trail in the non-stationary process, we linearize 
Eq. ( llOal for relatively small density disturbances, 5n — 
5N/Nq = n — 1. In accord with our theory and simulations, 
we also assume that the major gradients beyond the trail are 
perpendicular to the magnetic field, i.e., are directed mainly 
along the x coordinate. As a result, we obtain from Eq. i ll Oat 
a simpler equation, 



dtSr 



DdiJn = Ddt x ^ cs (x,t). 



(68) 



Solving this linear equation via the proper Green function and 
performing integration by parts, we obtain 



8n(x, t) 



1 



dt 



■ I '-I to (t-i} 3/2 

(x — x) 2 



x (x — x) exp 



AD(t - t) 



9 s </) res (i,t)£ii . 

(69) 



This expression shows that general density disturbances are 
determined by the entire distribution of the residual electric 
field over x from all previous times t. It is straightforward to 
check that for stationary 1OS and t to, Eq. ( I69l l reduces to 
dependence Sn(x, t) = — 4> rcs (x, t) corresponding to the local 
Boltzmann distribution for <f> ICS <C 1. 

Before applying a specific model for <9£<^ res (i, t) we note 
the following. When the local coordinate x varies from to 
00, the residual electric field — dz4> les (x, t) changes its sign. 
In the central region it draws plasma to the center, while out- 
side the central region it repels it. At any previous time t < t, 
the local density disturbances are determined by integral con- 
tributions from the two competing regions that generally do 
not balance each other. In the absence of the exponential fac- 
tor in the RHS of Eq. ( I69> . according to Eq. iD2\ . the two 
contributions would exactly compensate each other. The ex- 
ponential factor, however, breaks the balance. For example, in 
the center, x — 0, the attraction always dominates, so that the 
total density disturbance there is positive (plasma compres- 
sion). In the adjacent regions, the situation is opposite, so that 
the repulsion dominates there, resulting in plasma depletion. 

The strongest density disturbances are within the near zone, 
where the potential reaches its absolute minimum. For simple 
estimate of the potential (p in the early stage of trail diffu- 
sion, p ^ 1, we choose its zero-order approximation, (p(°\ 
Eq. (|56}, so that res w 0° cs w ip^ /2(1 + ip). Recalling the 
definition of the renormalized coordinate 77, Eq. d24bb . and 
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renormalizing the local time variable t, we obtain 



dx4>res (x, t) = — X 



where 



t 



77 if 

; (0 if 



> m 



4> 



l + V 



(70) 



(71) 



We restrict our analytic estimate to the center, x — 0, where 
we expect the largest positive density disturbances, and as- 
sume t S> to. As a result, performing the integration over 
space, we obtain from Eqs. d69l and d70l 



5n(0,t) 



1 



(7 1 -7 2 )dr, 



(72) 



where 



h = 



27 7 i 



6(l-r) s 



exp 



erf 



7/2 T 



4(l-r) 



(73a) 



2^T 



■ exp 1 



3(l-r) 



3/2 



4(l-r) 



(73b) 



are positive functions of r and I2 (r) < Ji (r). 

For small 771, which for 770 ~ 3-5 is possible only at suffi- 
ciently high altitudes where ip <C 77^ 2 , the two functions are 
mainly localized near r = 1 (t ~ t). As a result, we obtain 
Jq 1 Jidr ps and J Q T^dr sa rjl*Jir/?>, so that in this 

case <$n(0, t) ~ 77q/50. Under real physical conditions, the 
value of r/o can reach several units, so that the maximum den- 
sity disturbance of the background plasma in the near zone is 
of order unity. The largest value of density disturbances ob- 
served in our simulations was Sn(0, t) ~ 2, which is larger 
than those following from the linearized Eq. d68i but is of the 
same order of magnitude. 

At lower altitudes, if> > 1, according to Eq. (17 1 i . the value 
of 971 is not small but is always less than 770. As the value of 
771 > 1 grows, the maxima of I12 shift to smaller values of r. 
This means that the density disturbances at a given time t be- 
come more affected by electric fields at past time, t < t. In ad- 
dition, the ratio I 2 / 1\ becomes smaller, which means that the 
effect of the oppositely directed electric field beyond the ma- 
jor region of the potential distribution, lo/^/t > 771 , Eq. il 01 . 
becomes less important. For 771 ^> 1, J Q l\d,T ~ 4-^/7? In 771. 
The value of 771 < 770, however, can reach 3^4 as a maximum 
(at sufficiently large %\) corresponding to lower altitudes). In 
this case, the integral (Ji — I2) dr < Iidr can reach 
several units at most, so that Sn(0, t) should remain small. 

These simple estimates confirm our numerical observations 
that even in the strongest cases the relative plasma density 



disturbances in the near zone are of order unity at worst, see 
Fig- El In a H other situations or locations, such as for p > 1, 
t ~ to, etc., and especially in the far zone where the major 
current closure takes place, the density disturbances are much 
smaller. This shows that our basic theoretical scheme which 
relies on the Laplace equation for the electric potential in the 
far zone and neglects density disturbances is a reasonable ap- 
proximation. 



Weakest case, p> 1: adjustment of p for residual potential 



The simulations show that Eq. d63l > describes well the evo- 
lution of the simulated trail density for all parameters p, while 
the residual potential solution described in Sect. IV Bl shows a 
discrepancy for the late diffusion stage when p is large. To 
fit the evolution of the simulated potential, the expression for 
the master parameter p given by Eq. fl63i needs an adjustment 
p(t) — > pit) corresponding to an effective time lag. We argue 
that small background density disturbances beyond the trail 
are responsible for this discrepancy. 

To estimate background plasma density disturbances in the 
later stage of trail diffusion, we will use the zeroth-order ap- 
proximation for the residual potential given by Eq. d30i . In 
the limit of p ^> 1, when p(t) ~ yt, we have in the near zone 



where 



<^esOM) 



I{z) 



i(x/ypt) 

2(1+^)' 



(74) 



Substituting Eq. (174b for x = into Eq. ( I69> . we obtain 

1 

Sn(x, t) 



2(1 + ^)P(*) 



(75) 



Thus for p 3> 1 the background density disturbances are 
small, justifying the use of Eqs. (I68> and j69l >. 

We will use Eq. d75i to reconcile the discrepancy described 
above. One of the key factors in the derivation of the govern- 
ing equation for the near-zone potential, Eq. d23l >. has been 
the calculation of the density disturbance integral over the ef- 
fective ^-region in the near zone, AjVf, , as described in Ap- 
pendix|C] This effective region includes the trail with a nearly 
constant (p ICS and an adjacent region where a noticeable £- 
derivative of <p res builds up, as seen in Fig. [5] at y ~ —10. 
This derivative defines the fields and current closure structure 
in the far zone. The effective integral AN ( ^ 1 should include 
both the trail density and the disturbances of the background 
plasma. The calculation of Appendixlcltook into account only 
the former and completely neglected the latter. For a dense 
trail with p« 1, this proves to be justified because the rela- 
tive contribution of the background plasma disturbances to the 
integral proves to be small compared to the contribution of the 
trail density. On the contrary, for a much less dense trail with 
p ^> 1, the contribution of background density disturbances 
becomes comparable to that from the trail density. 
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In Appendix |F| we have obtained the effective integral 
AJV^j (^) and the relation between p and p in terms of two 
coupled parameters of order unity, j3(p) and $(p), Eqs. ( IF5> 
and JF6i . Being unable to obtain these parameters analyti- 
cally, we have used our simulations for tp > 0.05 to obtain 
approximate Eq. JF8i . For smaller tp, we proposed Eq. iF9l . 
which represents a conjecture and needs a special study. 



or, in the explicit form, 



Pit) 



W ( 1 /, e ( 1+ ^^*+'/') - tp 
l + ip 



(79) 



where W(x) is the Lambert W-function. 

To calculate the residual potential for sufficiently large p, 
the parameter p needs an adjustment, p — > p < p, 



VI. SUMMARY OF ANALYTICAL RESULTS 

In this Section, we summarize our major analytical results 
which can be directly applied in comparisons with observa- 
tions. To use the theoretical results for practical applications, 
we need to pass from the renormalized variables of Sect.lVlto 
the original ones. Where appropriate, we will use simplified 
versions of analytic expressions. 

The original residual potential <p KS (x,y,t) is defined in 
Eq. @. According to Eq. J24i . the original residual poten- 
tial in the near zone, (^(x, t) = <p Tes (x, 0, t), in terms of the 
actual coordinate x and time t, is given by 



i 



7 



x, p(t) 



(76) 



2(1 + V) r VV Dp(t) 

where 7 is defined by Eq. ( I26l l. Assuming the E-region con- 
ditions, v en ~ lQvin and m, ~ 30m p , we write 7 in a 
form convenient for comparison with our FlexPDE simula- 
tions, Sect. II VI where n = Nxra,a(to)/No and D = 1, 



1 



37.1 



7 



(77a) 



In terms of the trail line density Nu n , background plasma den- 
sity No, temperatures T e ^, and geomagnetic field Bq, we ex- 
press 7 and the diffusion coefficient D as 



7 = 



2ttNq {T e + Tj) 
^ (1 + V) N lin eB 

1.08 xl0~ 2 (T e +Ti\ /0.5xlO- 4 T 



vAKi + VO V ioook 

iV \flQ 14 wT 1 
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lin 
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T e + T t 23.2 / T e + T t 



rriiVin 



^P V 1000K 



0.5 x 10~ 4 T 



(77b) 



2 — 1 
m s . 



(77c) 

The only time dependence in </>° s is associated with the di- 
mensionless parameter p(t). The function <p(f),p) has dif- 
ferent approximate expressions depending upon the range of 
p, as described in Sect. IV Bl This parameter monotonically 
varies from small values in the early diffusion stage to large 
ones in the later diffusion stage. According to Eqs. J62I and 
A63I . the time dependence of p can be approximately deter- 
mined by 



it = P + 



1 



1 + tp 



In 1 



i + y> 

tp 



(78) 



where 

hp) 



p(p) = §(1-/5(p) + V /i -2/3(p) 



0.92p 



Hp) 



(l + V)(2.2 + p) V TT/O 

p fjt 

2(2.2 + p) V p 



(80) 



— , for V> 0.05, (81a) 



for ip< 0.05, (81b) 



and the relationship between jt and /? is given by Eq. (178b . 
Equation ( 181 al l includes the approximation for /3(p), Eq. 1 IF8I 1. 
obtained by fitting FlexPDE runs for t/j > 0.05, while 
Eq. J8 lbb is an extrapolation to smaller tp. 

In the earlier stage of dense-trail diffusion when p(t) « 
V>7*/ (1 + ip) <C 1, a simple approximation for the near-zone 
potential is given by the piece-wise function, Eq. ( I36> . 



1 



2(1 + V) 



Zk + IL. if ^(^t)! <7?0 , 

^ if |7 7 (x,t)| >77 , 



(82) 

where 77(2;, t) = x^/D/jft)] 1 / 2 and?/ - 2 [ln(l/p(t))] 1/2 . 

For simple estimates of the residual potential near the trail 
at a later stage, p(t) ~ jt > 1, one can use the zero- 
order Eq. d30i . To compare with simulations, however, we 
should use the more complicated, but more accurate, first- 
order Eq. J3 II . and replace p by p, Eq. d80i . and r\ by fj(x, t) — 
xlj/Dp(t)} 1/2 , 



1 



2(1 + V)p(t) 



S(r)) 



V / l+4A(/5) 



4A(p) 
l + 4A(p) 



(3-2A(/5)7) 2 )5(7- y Vl + 4A(/5))]}, (83) 



where the functions 5, A, and A are defined by Eqs. JC 17i and 
(I34> . Given 4>® es (x, t), the residual potential in all locations is 
given by Eq. ( IB4> . At large distances from the trail, x 2 + 
a 2 y 2 3> Dt, according to Eqs. dB5b and d25l >. the residual 
potential has a 2D quadrupole structure, 

CiipjpDj^-aY) 
^ V > t] M 2(l + ^)7(^+aV) 2 ' (8 } 

where Ci(p) = lim^|_ >00 [r7 2 ^(77,p)], Eq. {D3}. 
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The expression for the trail density is much simpler than 
those for the potential. According to Eq. ( 149 \ with q = p/yt, 
Eq. \25\ . for all values of p it has a nearly Gaussian form, 



AriTrailO^y,*) ~ An 0l 



' U)P0_ 

tp(t) 



exp 



1 fix 2 y 2 
AD \p{t) + t 



(85) 

where p = p(t ) w 7^0/(1 + t/>), Eq. d58l . and no adjust- 
ment for all p. 

Eliminating t from Eqs. J78> and (I85> . we express the peak 
Gaussian density, Anp ea k = AnTr a ii(0, 0, t), in terms of p(t), 



An Pcak = < 2604/ p 



(1 + V) p + In ( 1 + p 



(86) 

while the same quantity in the self-similar solution is given by 



An ss 



Peak 



20o 



P' 



1 infi + i^p 



1+V 



(87) 

These expressions have a universal form independent of the 
initial conditions, provided t > 2t and An(t ) ^> 1. 



VII. COMPARISON OF THEORY AND SIMULATIONS 

In this Section, we compare our analytical theory with sim- 
ulations outlined in Sect. IIVI We start with comparison of 
our simulations with the theoretical expression for the trail 
density, Eq. ( 1851 . In the early stage of dense-trail diffusion, 
p -C 1, the nearly Gaussian peak in simulations closely fol- 
lows the self-similar solution (SSS), Eq. (I47> . though it shows 
a slightly faster diffusion, as seen in Fig. UOt a). The density 
predicted by Eq. (1851 is closer to the SSS curve than to the 
numerical one because the interpolation expression for p{t) 
given by Eq. (I63> is less accurate for p <C 1 than it is for 
P>1- 

When p ~ 1, the analytic solution for the density peak 
given by Eq. ( I85> starts deviating from the SSS and becomes 
closer to the numerical solution, as illustrated in Fig. UOt b). 
When p becomes larger, p ^> 1, the analytical theory shows 
an excellent agreement with simulations, while the SSS pre- 
dicts a noticeably slower diffusion, as shown in Figs. UOf c) 
and (d). 

Figure ^2 shows the peak trail densities vs. p taken from 
the simulations, analytical theory, and self-similar solution. 
The numerical solution is shown by separate groups of points 
taken from several different runs (each group has its own point 
shape). Each run started at to = 1 with different values of the 
initial peak trail density. In each group, consecutive points 
correspond to equidistant moments of time: t = 1,2,3,.... 
The first point of the group always lies on the dashed curve 
corresponding to Eq. ( 1871 because the SSS was the initial con- 
dition for each run. However, starting from t = 2, the nu- 
merical points closely approach the theoretical curve given 



by Eq. d86l >. while the SSS solution given by Eq. d87i re- 
mains noticeably offset. The theoretical curve overlays the 
numerical points for p > 1, while for p <gc 1 it shows a 
slight deviation from the numerical points (see the beginning 
of this Section). The transition from an anisotropic diffu- 
sion to a more isotropic one occurs near the inflection point 
about p = 1. Notice that this takes place when the trail peak 
density remains well above the background plasma density, 
An max ~ (OaViiy 1 ^ 80 for $ = 0.05. 

Now we compare with simulations the theoretical expres- 
sions for the residual potential in the near zone. Figure [T2l 
shows the residual potential along x for the same conditions 
as in Fig.no| Because the parameter p spans a broad range of 
values from small to large ones we will apply either Eq. fl82l > 
or Eq. ( 1831 . Figure [21 a) shows that for p -C 1 the simple 
piece-wise approximation agrees reasonably with simulations 
in all areas not too close to the two positive bumps of the po- 
tential. Indeed, while there is a significant difference between 
the values of the potential minimum at x = 0, the potential 
derivative (the residual electric field) is the same in the inner 
region characterized by the parabolic dependence and occu- 
pied by the trail, see Fig. UOf a). On the other hand, a good 
agreement also exists well beyond the trail, where the resid- 
ual potential decreases with increasing |x| and the correspond- 
ing electric field changes its sign. The zero-order piece-wise 
approximation is rough in the transitional zone near the two 
potential maxima, where it has a discontinuity in the electric 
field. The maximum electric fields in the piece-wise formula 
are reached near the discontinuity points, approaching them 
from inside. The maximum electric fields in simulations are 
reached at some locations in the inner region closer to the cen- 
ter and hence have smaller values. Thus the simple analytical 
formula yields nearly correct electric fields everywhere except 
the transitional zones between the inner and outer regions, 
where it overestimates the electric field magnitude. We have 
attempted modeling the transitional electric field with higher- 
order interpolations to provide a smooth transition, but this 
underestimated the field. The error for the maximum electric 
field, however, remains within the range of tens percent for all 
our simulations. 

For p ~ 1, the theoretical expression given by Eq. (1831 with 
unadjusted p(p — > p), agrees well with simulations practically 
in all locations, as seen in Fig. I12f b). As p becomes large, 
Eq. (1831 with unadjusted p shows a significant discrepancy, as 
seen in Fig. 1121 c) and (d). However, if we adjust the parameter 
p, p — > p based on matching of the potential minima for the 
two solutions and apply Eq. J83i then the discrepancy prac- 
tically disappears. It is important that the analytical solution 
with only one adjusted parameter p — ► p causes the theoret- 
ical (j>® cs to match the numerical solution well not only near 
the potential minimum but everywhere. The relation between 
p and p based on our simulations results in the empirical p- 
dependent adjustment coefficient j3, Eqs. d80t and J8 U . We 
emphasize that the adjusted parameter p is only needed for 
the residual potential and not for the trail density, as described 
in Sect. lVTjl 

In this Section, we showed the comparison of theory and 
simulations mainly for small ip. For large values of ip, the trail 




FIG. 10: Trail density distribution along x for ip = 0.2. (a): An = 10 4 , t = 10, p « 1.4 x 10~ 3 . (b): An = 1000, t = 20, p » 0.5. (c): 
Ano = 100, t = 10, p ~ 5.37. (d): Ano = 10, t = 20, p ~ 160. Solid curves: the numerical solution; dot-dashed curves: the self-similar 
solution, Eq. <47t : dashed curves: the analytical solution, Eq. 1851 . 




diffusion in all stages of trail diffusion was nearly isotropic 
with small values of the residual electric field, in full accord 
with the analytical theory. 



VIII. DISCUSSION 

In this section, we will start by discussing some caveats, 
then we will estimate the induction electric field in the plasma 
trail, and finally we will dwell on plasma instabilities. 



A. Caveats 



FIG. 11: Peak trail density vs. parameter p for tjj — 0.05. 
Points shown as alternating crosses, diamonds, boxes, and cir- 
cles represent the numerical solution from several runs: Ano = 
10 4 , 10 3 , 100, 50, 30, 10. For each run, the consecutive points (from 
left to right) correspond to t/to = 1, 2, 3 . . .. Solid curve shows 
the analytical solution given by Eq. 1861 . Dashed curve shows the 
self-similar solution given by Eq. 1871 . 



In our theory, we have made a number of assumptions 
which are not perfectly valid. Among those were assump- 
tions about constant electron and ion temperatures and about 
one sort of ions. However, the initially hot temperatures of the 
newly produced meteor trail plasma need some time for cool- 
ing. Also, this plasma includes material different from the 
ambient atmosphere, so that there may be at least two kinds 
of ions with different masses. We will include these factors in 
future work. 

We also note that our full analytical theory has been de- 
veloped only for the particular case of a mutually orthogonal 



Residual potential Residual potential 




FIG. 12: Residual potentials at y = 0, 4>® cs (x, t), corresponding to density distributions in Fig. 1101 Solid curves show the numerical solutions, 
while dashed and dot-dashed curves show the analytical solutions according to different equations, (a): Ano = 10 4 , t = 10, p ~ 1.4 x 10 -3 , 
dashed curve corresponds to Eq. 1821 . (b): Ano = 1000, t = 20, p ~ 0.5, dashed curve corresponds to Eq. 1831 with p — » p. (c): Ano = 100, 
t = 10, p as 5.37, dashed curve corresponds to Eq. I83t with p — > p, dot-dashed curve corresponds to Eq. 1831 with p ~ 3. (d): Ano = 10, 
t = 20, p » 160, dashed curve corresponds to Eq. <83l with p — > p, dot-dashed curve corresponds to Eq. <83l with p ~ 68. 



meteor trail axis and magnetic field, 9 = 90°. A more general 
situation occurs when the angles between the trail axis and 
Bo, 6, are arbitrary but satisfy restrictions given by Eq. Q. 
It is only discussed for the self-similar solution described in 
Appendix 1X1 If we ignore the effect of electron Hall currents 
into the process of the meteor trail diffusion then we can apply 
all our results where according to Eq. (1121 the quantity a~ x , 
Eq. ( IB2I . and the numerator of p, Eq. ( I25> . acquire an addi- 
tional factor sin 9. Our preliminary numerical computations 
show that the electron Hall current at 6 ^ 90° affects the me- 
teor trail diffusion in such a way that the trail density contours 
become more isotropic and rotated at a small angle in the x, y 
plane, as the SSS does (see AppendixIXl. 



B. Induction electric field 

Now we estimate the induction electric field associated with 
the meteor plasma trail currents and show that, except for 
extraordinarily high-density meteors, this field is negligible. 
This is of importance for the electrostatic field approxima- 
tion employed in our approach. The induction electric field 
may only occur due to disturbances of the magnetic field, <5B, 



which in turn are due to the electric current, j, caused by the 
plasma trail diffusion. Note that the magnetic field distur- 
bances, even in spite of their extreme weakness, can be ob- 
servable using a sophisticated measurement technique for di- 
agnostics purposes. Furthermore, the initial formation of the 
trail current may have caused ELF/VLF signals observed on 
the ground and correlated with the meteor showers II241 12511. 
According to the Maxwell equations, we have 

V x E = - d t SB, (88a) 
Vxffi-M (88b) 

where fio is the permeability of free space and we neglected 
the displacement current (all characteristic speeds in the trail 
diffusion are much less than the speed of light, etc.). The net 
electric current is due to the fact that electrons and ions have 
different responses to the external force. Being proportional 
to the plasma density, the net electric current is mainly con- 
centrated within the trail. It is predominantly formed by the 
unbalanced E x B drift of electrons. In the strongest case, the 
trail density behaves roughly in accord with the self-similar 
solution, so that for a simple estimate we can use the corre- 
sponding expressions for the current given by Eq. (IA10> . This 
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current with the density^ = (eSl e xs\Ti9 /2v en t)N &s , N ss = 
(C*/i)exp{- [(1 + xj))x 2 /ijj + y 2 ] /ADt), Eqs. JXTJ and 
(IA6> . is directed along the trail axis and flows in opposite di- 
rections in the two halves of the trail, x > and x < (in the 
actual 3D, spatially inhomogeneous, and restricted trail, the 
current forms a closed loop). Setting as above 9 = 90° and 

1 /2 

taking a typical scale along x as Ax = [A.ipDt/ (1 + ip)] 
and N ss ~ AT max , according to Ampere's law, Eq. (I88b> . we 
estimate the typical magnetic field disturbance as 

2ep n e D { ip \ M 2/i (T e + r t ) 

d# — 7V max = — — — A max , 



1 + 



(1 + So 



(89) 

where in the last equality we used the definitions of D, 
Eq. il It . and ^, Eq. (|2j. The relative magnetic field distur- 
bance is 



SB p N n 



(90) 



where /3q is the ratio of the undisturbed total plasma pressure, 
Nq(T c + Tj), to the magnetic pressure, Bq/2/xq, 



Po « 1.39x10" 



An 



10 n m -3 M 10 Q0K 



T e + Ti\f 0.5 x 10~ 4 T 



S 



In the E-region ionosphere, ATq ~ (10 9 -10 11 ) m 3 , T e ~ 
Ti ~ 300K, B ~ (0.25-0.6) x 10~ 4 T, so that the relative 
disturbance of the magnetic field is small, unless the plasma 
trail is extremely dense, N max /N > 10 6 -10 8 . 

To estimate the contribution of the induction component 
into the total electric field, we can estimate the ratio of | V x E | 
to |V • E| > (T e + Ti) /e{Ax) 2 . According to Eqs. J88H and 
( I89> . after simple algebra we obtain 



V x E 



V • E 



< 



46n 



SB_ 

(1 + V) Z V-Bo 



(91) 



This ratio has an a additional factor in front of SB /Bo, which 
is always small since according to Eq. 0, 40oVW(l + 
■0) 2 < 0.3256 ^ 0.44 x 10~ 2 . Because SB/B in the reg- 
ular meteor trail is small, the induction electric field proves to 
be even smaller. This justifies our initial assumption that the 
total electric field is nearly curl- free, E = —~V<p. 



C. Dynamics of electric field and plasma instabilities 

The ambipolar electric field associated with trail diffusion 
may drive plasma instabilities responsible for observable non- 
specular radar echoes. In this Section, we will make simple 
estimates of the driving field and instability threshold. 

The driving field is determined by the total external force 
acting on electrons, E rcs = — (T e + T;)V</> rcs , where <fi Tes is 
the residual potential defined by Eq. Equations J8 21 to 
J84i give approximate analytical expressions for the residual 
potential if the magnetic field and the meteor trail axis are mu- 
tually orthogonal. These expressions depend upon the dimen- 
sionless parameter p(t) given by Eq. M9\ in terms of the ef- 
fective rate 7, Eq. \2bl . The latter parameter is proportional to 



an effective trail-background interaction cross-section a c s = 
N\i n /No- After the critical time given by Eq. j66l >. the diffu- 
sion process becomes more isotropic and the residual poten- 
tial decreases drastically. The critical time i cr varies depend- 
ing upon the ionospheric conditions and meteor parameters. 
Nighttime conditions with low Ao and given 2Vu n are equiva- 
lent to daytime conditions with much larger No and the pro- 
portionally increased column line density Ny m . For example, 
assuming equatorial day-time ionosphere, No ~ lO 11 !]! -3 , 
and a typical linear trail density, Au n ~ 10 14 m _1 (J, 0], 
we obtain a critical time of tens of milliseconds. For mid- 
or high latitudes, night-time conditions, No ~ 10 9 m~ 3 , and 
Min 10 15 m _1 (or faster and/or bigger meteoroids), we ob- 
tain that the critical time may reach tens of seconds. 

The strongest electric field occurs in the early diffusion 
stage, p -c 1, and at the edge of the nearly parabolic re- 
gion of the potential, Eq. fl35t . Using Eqs. A36I . j46> . and 
Eq. J24> . we obtain the maximum value of the residual electric 
field, |S™ ax | ~ [D~fln(l/p)/p] 1 / 2 m i v in /e(l + tp), where 
the diffusion coefficient D is defined by Eq. ( II It . In the 
later stage of trail diffusion when p ^S> 1, Eq. ( I27l > gives 
[d v <p(T])] max — 0.7/ p, ignoring the adjustment of p described 
in Sect. IV Dl In the original variables, we obtain |£"™ ax | — 
0.35(Dj/p) 1 ^ 2 rriih'i n /e(l+il))p. Interpolating between these 
limiting expressions, we can write a simple formula, 



e(l + V) 



D7 
~2p~ 



In 1 



0.3 

_2 



1/2 



(92) 



roughly valid in the entire domain of p. 

Now we estimate the Farley-Buneman (FB) instability cri- 
terion For a homogeneous plasma, the simplest 
FB instability criterion, obtained using the two-fluid plasma 
model for sufficiently long-wavelength waves, is given by 
I Vo| > (1 + ^)C S El, where V = E rcs x B /B 2 is 
the E rcs x B drift velocity and C s = [(T e + T^/rrn] 1 ' 2 = 
{Duin) 1 / 2 is the ion-acoustic speed. Applying this criterion 
to the maximum field given by Eq. ( I92K expressing the cor- 
responding drift speed as |Vo| = e|-E™ ax |/mi$!i, and us- 
ing Eqs. (|2j and 0, we write the FB instability criterion as 
|£™ ax | > (1 + ^)( J Di/ in ) 1 / 2 m i n i /e. Expressing E™ 3 * in 
A^in, we reduce the FB instability criterion to 



0.3 
"^2" 



In ( 1 + ^ ] > P, 



(93) 



where 



7r(T e + T l )m i No 



1000K \ / S, 



T e +TWV0.5 x 10~ 4 T J \10 14 m 



Min \/10 n m- 3 



No 



(94) 



If this criterion is satisfied when the trail initially forms, then 
the instability starts generating plasma irregularities. If the in- 
stability persists for the sufficient time, then turbulence will 
develop and partially saturate through nonlinear processes. 
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FIG. 13: Altitudinal dependence of FB instability duration £fb for the equatorial E region during (a) daytime (No = 10 11 m 3 ), where curve 



1 is for Nuu = 10 m and curve 2 is for Nim = 10 m , and (b) nighttime (No 



ioV- 



), where curve 3 is for Min = 10 m and 



curve 4 is for Nu n — 10 J m . Curves 1, 2, 3, and 4 correspond to <j c ff = 10 , 10 , 10' , and 10 m , respectively. 



Because pit) monotonically increases with time, see Eq. d63l >. 
then at some moment, t = tpB, the two sides of Eq. d93i 
become equal. At this moment, the linear growth of the FB 
instability starts being replaced by linear damping and irregu- 
larities will diffuse away. At t < £fb, the linear FB instability 
sustain plasma turbulence at a certain level, while at t > £fb 
there is no more free energy to sustain the turbulence, so that 
irregularities will quickly disappear. 

The instability duration £fb depends critically upon the al- 
titudinal parameter ip and the effective trail-background inter- 
action cross-section a e g = Nu n /No. Due to this, the night- 
time conditions with low No will produce longer lived meteor 
trail than will the daytime conditions. Figure IT3l shows the 
altitude dependence of the instability duration tpB for several 
constant <j c & during daytime and nighttime conditions. Notice 
clear peaks of t-pB at some intermediate altitudes which in- 
crease with a c g. The non-specular echo boundary for a given 
trail, like that in Fig.^ should roughly follow the altitude de- 
pendence of tpB- However, because column plasma density 
varies along the meteor trail and due to other inhomogeneities, 
we expect more variability than this model predicts. Mea- 
surements of the evolution of the trail echoes, in combination 
with other observations, should enable us to retrieve useful 
information about meteors, ionosphere and atmosphere. Im- 
plementation of this procedure requires better models of in- 
stability generation than that used above, as well as models 
of ablation and ionization to give improved estimates of Nu a 
id. 



IX. SUMMARY 

In this paper, we have described analytical theory and finite- 
element simulations of trail diffusion and fields for the mutu- 
ally orthogonal trail axis and magnetic field. Unlike previous 
models, this theory includes both the trail and the background 
ionospheric plasma. This has two major effects: (1) a natural 



restriction on the ambipolar electric field that otherwise would 
infinitely grow with distance from the trail and (2) a later-stage 
transition from sharply anisotropic (for tp <C 1) diffusion to 
nearly isotropic one. The former is important for plasma in- 
stabilities responsible for non-specular radar echoes, while the 
latter is important for interpreting specular radar echoes. 

A key element of the present treatment is the introduction 
of the residual potential, <fi rcs , defined by Eq. 0. Its gradient 
describes the total force acting on electrons. Due to high elec- 
tron mobility along the magnetic field Bo, the typical scale of 
4> Tes spatial variations in this directions are much larger than 
that in the perpendicular direction, while the typical scales of 
trail density variations in both directions are comparable. This 
fact, which had not been realized in earlier simulations of me- 
teor trail diffusion, requires setting the simulation box bound- 
ary along Bo far from the trail boundary and simultaneously 
resolving the two different scales parallel to Bo- In our sim- 
ulations, we have overcome the computational difficulties by 
employing a finite-element software FlexPDE with the adap- 
tive cell structure. 

Based on the insight from simulations and using the large 
ratio of the electron and ion mobilities along the magnetic 
field, we have developed an approximate analytical approach. 
In this approach, the problems of trail diffusion and of spa- 
tial distribution of the residual potential are treated separately, 
while the coupling between the two is made via parameters 
and approximate solutions. In particular, Gaussian approxi- 
mation of the trail peak has allowed us to to treat the resid- 
ual potential in terms of one parameter p. This parameter 
is proportional to the square of the trail peak dispersion, as 
described by Eq. j25L and monotonically grows with time. 
Due to high electron mobility along the magnetic field, the to- 
tal force acting on electrons in this direction is much smaller 
than the corresponding components in other directions. For 
the particular case of mutually orthogonal trail axis and the 
magnetic field, this has allowed us to reduce the original 2D 
description to a ID linear integrodifferential equation given 
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in two different forms by Eq. (1231 . This governing equation 
is for a residual potential in the near-trail zone and it has a 
unique solution depending upon p. The approximate solu- 
tion of Eq. d23i depends on the range of p, as described in 
Sect. IV Bl Using this solution, we have obtained the approx- 
imate expression for p(t), Eq. d63i . which closes the entire 
description of trail diffusion and fields. Note that the near- 
trail potential, which is easily spread along the magnetic field, 
may create significant disturbances of the background plasma 
beyond the trail, as described in Sect. IV Dl 

Comparison of analytical theory with simulations have 
demonstrated good agreement between the results with one 
exception. At a later stage of trail diffusion, the parameter p 
for the residual potential should be replaced by an adjusted pa- 
rameter p as described in Sect. lVIll We have identified the na- 
ture of this deviation analytically in terms of the disturbances 
of the background plasma, but to quantitatively relate p and p 
we have invoked simulations. 

The analytical theory and simulations have allowed us to 
estimate the spatial distribution of the ambipolar electric field 
within and near the plasma trail. These fields are crucial for 
plasma instabilities responsible for generation of field-aligned 
electron density irregularities observed by high-power large- 
aperture (HPLA) radars as non-specular echoes. Measuring 
the characteristics of non-specular echoes and some other 
characteristics should allow one to retrieve an important in- 
formation on the meteoroids and the surrounding atmosphere. 
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diffusion starts from an infinitely thin and dense plasma col- 
umn with a given line density. The electron Hall velocities 
give rise to the non-diagonal coefficients, A xy = B xy . In- 
equalities given by Eq. ( IA3I mean that the contours of the 
constant density and the residual potential form ellipses in the 
xy-plane, whose major axes are rotated with respect to the x 
and y axes through a common angle x determined by 
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The constant C in ( IAU is expressed in terms of the conserved 
linear density (along the z coordinate), 2Vu n , as 
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For arbitrary electron and ion mobilities it is hard to obtain ex- 
plicit analytical expressions for the coefficients Aij and £>y . 
However, under restrictions described by Eq. H2i . to lead- 
ing order accuracy with respect to the small parameters 
(ipQ)" 1 , we obtain 



A xx — 1 + B x 



B xx — 



sin 



A 



yy 



1 + Byy , Byy 



4> ' 
i + cos 2 e 



e 2 



1 + cos 2 6 



(A6) 



A — R 

s*-xy — ^xy 



2/i sin 2 6> 



2fii 



• cos 6. 



The non-diagonal coefficient, A xy = B xy , is always small 
compared to A xx , and A yy w 1, but not necessarily to 
B xx and B yy <C 1. The rotation angle, Eq. dA4l >. is small 

| tan 2%| « 2u en cos 6/(fl e sin 2 9). Because 



APPENDIX A: SELF-SIMILAR SOLUTION 

In this Appendix, we obtain explicit expressions for the 
self-similar solution (SSS) proposed in the general form, but 
not found explicitly, by Jones 1 17]. In our notations, this so- 
lution (denoted below by superscript 'ss') follows the ansatz 

n(r, t) = n ss (x, y, t) and ros (r, t) = (j>f es (x, y, t), where 
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n ss (x,y,t) = — exp 
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C s (x, y, t) = n, -' :1 ' ' + const, (A2) 

with positive diagonal coefficients An, Bn and 
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This is the solution to Eq. (II 01 in an infinite and homogeneous 
neutral atmosphere with no background plasma, provided the 



The residual potential is stretched along the coordinate y in 
accord with the qualitative discussion in Sect.||l] which holds 
under condition 



(A8) 



This restriction due to the electron Hall velocity is stronger 
than that of Eq. ( fTal . For the trail strictly perpendicular to B, 
cos 9 = 0, p = 0, Hall velocity is directed along z and plays 
no role in 2D trail diffusion. 

To the same accuracy, the particle fluxes are given by 
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T iz = 0, T e2 = - E^^i y> t ) ( A9) 

"en ^ 

In the SSS, the flux components in both x and y directions 
are equal for electrons and ions. The only disparity is in the 
flux component along the z-axis due to electron Hall velocity. 
The net electric current is directed along the trail axis with the 
current density 



Jz 



x sin ( 
2t 



n ss (x,y,t). 



(A10) 



APPENDIX B: RESIDUAL POTENTIAL IN THE FAR ZONE 

In this Appendix, we solve for (p rea in the far zone |£| 3> 1, 
where £ = yf(Dt) 1 / 2 is a coordinate parallel to Bo- This 
zone is located well beyond the plasma trail. Neglecting 
plasma density disturbances allows us to reduce Eq. dl5b> to 
the Laplace equation 

<%0rcs + Secures = 0, (Bl) 

where we have renormalized the variable £ as £' = ct£ with 

, , \ 1/2 

_< ^ =e (l + ^) 1/2 «l. (B2) 

In terms of £ and £', the near zone described in the following 
Appendix, reduces to a thin strip extended along £. Assuming 
that potential disturbances have typical scales of order unity 
with respect to both coordinates, we can approximate this strip 
by a cut in the £. £'-plane at £' = ±0. Here the signs ± mean 
infinitesimal offsets from zero to either positive or negative 
directions. Because of the mirror symmetry with respect to 
the £-axis, Eqs. dl9b> . we will consider only the positive half- 
space, £' > 0. 

The nearly constant value of the residual potential in the 
near zone represents a boundary condition for ra s(£, £') at 
the cut, 



fe»&C)] C '=±o«<.(0. 



(B3) 



Other boundary conditions are given by </>res(£>C') ~~ * as 
£, £' — > ±oo, The solution of Eq. (IB 1 i with these bound- 
ary conditions (the Dirichlet problem for the upper half-space, 
£' > 0) is given by 



^res(£,£') - 



c r c(T)dr 

* loo (r-0 2 +C 2 

, (r-£) 2 +a 2 £ 2 



(B4) 



In particular, for £ 2 + £' 2 > 1, Eq. (fB4l reduces to a 2D 
quadrupole, 



</>re S (£,C) 



(e + c 2 ) 2 ' 



(B5) 



In polar coordinates, r and 9 are defined as £ = rcos9, 
£' = r sin 0, the quadrupole potential is (j) Ies (r, 9) « 
(Ccos26')/r 2 , while the electric field lines of force are de- 
termined by (sin 29) jr 2 — const. Here the constant C is 
determined by the distribution of </> ros (£) along the cut £ = 0, 
as discussed in AppendixlDl 

Now we obtain general relations which follow from 
Laplace Eq. JB 1 i and will be used in Appendix Icl to derive a 
closed equation for </>r CS (£) when combined with the residual 
potential in the near zone. Introducing a complex coordinate, 



Z £ • *C, 



(B6) 



we consider £' > as a complex half -plane and introduce a 
complex potential, 



(B7) 



vP(Z) = rcs (£,£')+*^ (£,£'), 

where the function 14 7 (£, £) also satisfies Laplace Eq. JB II . 
This is an analytic function of the complex coordinate Z with 
</> res and W related by the Cauchy-Riemann equations, 



(B8) 



At the cut, £' = +0, we denote tu°(£) = W(£, +0), so that 
from Eq. (IB8i we have 

= dcW\ c=+Q , d^ ICS \ (=+Q = - d e w°. (B9) 

According to the mirror symmetry, Eq. dl9bl . the function 
0res(O i s even, while u>°(£) is odd, 

<PL(0 = <t>L(-0, *A0 = • (bio) 

At large Z, the residual potential <p ICS oc 1/\Z\ 2 decreases 
faster than 1/\Z\ as \Z\ — * oo. According to the Cauchy- 
Goursat theorem, § c i 5>(Z)dZ = 0, so that we have 



w°(£K = 0. 



(Bll) 



While the integral relation for u>°(£) is trivial because w° is 
odd, the integral relation for the even function 0r CS (£) repre- 
sents an important constraint. 

Now we proceed with the complex potential ^(Z). In 
the upper half-plane, we consider the continuous contour 
C, which includes the axis £' = +0 with an infinitesimal 
half-circle around Z' = £ + iO, and the infinite half-circle, 
\Z\ — > oo. Because ^ is an analytic function and the pole 
Z = Z' is beyond the area closed by this contour, the residue 
theorem yields § c ${Z)/(Z - Z') dZ = 0, or 



*(t) 
r-£ 



dr = iir^ (£) , 



(B12) 



where V denotes the principal value of the integral along the 
real axis. Separating in Eq. ( IB 121 the real and imaginary part, 
we obtain 



i r 

t>L(0 = --P 

7T J-< 



w°(0 = --V 

7T 



w° (r) 



d.T, 



dr. 



(B13a) 



(B13b) 
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These equations are equivalent to the well-known Kramers- 
Kronig dispersion relations in optics, plasma physics, etc., 
which are derived in the same way. Applied to an arbitrary 
function 4>® cs (£,) [or w° (£ )"], Eq. JB 1 3i is also known as the 
Hilbert transform 1 29, 30]. Equation (IB 13l > will be applied in 
AppendixOand Sect. lVBl 



APPENDIX C: NEAR-ZONE POTENTIAL 

In this Appendix, using the results of Appendix 151 we de- 
rive the governing equation for the residual potential in the 
near zone, |C| <C 1/a, where a = 0o(l + ip) 1 ^ 2 <§; 1 is de- 
fined by Eq. (IB2I . To obtain a closed equation for </> res (£), we 
write Eq. dl5> for /i = in a conservative form as 

dt (tAn) - l«N- 7 ^ 



V • [(1 + An) V</w] = 0, 



(Cla) 



d c [(1 + An)% 



1 

a 2 
V 2 An 
1 + 



= 0, 



(Clb) 



where An = AN/Nq is the total density disturbance which 
includes both the plasma trail and the background plasma dis- 
turbance. We assume that the major density disturbances are 
concentrated within the near zone, ( <C where the resid- 
ual potential is only weakly ^-dependent. In this zone, which 
overlaps with the far zone £ ^> 1, we represent the residual 
potential as 

4>r es (Z, C) = riZit + 0, M « \4>lZ\, (C2) 

where c/>rcs(? , C) is the component in the far zone and therefore 
satisfies the Laplace equation described in Appendix 151 Be- 
cause most of the potential changes occur in the background 
plasma (far zone), our assumption \6(f)\ <C l^resl is well jus- 
tified as can be seen in the example solution shown in Fig. [5] 
The function 0™£ (£ , () varies along £ with a large scale-length 
typical for the far zone, so that within the near-zone it varies 
approximately linearly 



^(e,C)«Cs(^o) + C9 c C I s (0, 



y rcs 



far , 



(C3) 



where |C<9c^'esl < mm(£, 0)|, with the derivative <9 C 0™' S re- 
maining nearly (^-independent across a transitional zone be- 
tween the near and far zones as seen in Fig.|5] The perturba- 
tion 8<j)(£, £), associated with strong density disturbances in 
the near zone, is relatively small, but its ^-derivative is not 
small compared to c^cA^s- Furthermore, it is this perturbation 
that provides a smooth transition from the zero ( derivative of 



</> res (£, C) at m e symmetry center, d^(p 



far 
res i 



= o, 



(C4) 



to the finite derivatives d^<p ICS « defiles in the far zone 
where 6<p(£, () gradually disappears. Now we subtract from 



Eq. (IClb> the Laplace equation for </>*cs an d integrate over £ 
from to £ = (i, where d -C a -1 is located in the transi- 
tional zone where both 5(f> and d^S<fr are zero. Also, at £ = £i 
we neglect the density disturbances An and the perturbations 
S(f) Tes with their derivatives, except in the term cx 1/a 2 . We 
will also ignore within the range |£| < d all density dis- 
turbances other than the trail plasma itself. (This introduces 
some error which we correct in Sect. IV D2l and Appendix|F|) 
The range |£| < d includes practically the entire trail, so that 
we can extend the upper limit of all density integrals to infin- 
ity. As a result, using Eq. ( IC4t . we obtain 



d s (AN c d ( c 

Here 

AiV fl (0 = 



> far l 

res ) 



far 
res 



AndC,^ AiVTrail(C) 



4 (AN ( ) = 0. 

(C5) 

An Trai i dC,. 

(C6) 

where Anxraii is the part of the total density disturbance as- 
sociated only with the trail plasma. 

To good accuracy, the trail plasma distribution is described 
by a Gaussian distribution, see Sect.llVland IV CI 

A (C /" -A n ° to I" ( ¥ . ^ 

An Trail( e, c,*) = — ] /^r ) ^[-{-^r ) + T 

(C7) 

where qo = ip/(l + ijj) and q(t) gradually varies between 
qo and 1. Because 0^ is the far-zone potential, we apply 
Eq. JB9b . d(<fi T % = — ac^ui ^). Bearing in mind Eq. dC3l >. 
in Eq. dC5l > we can eliminate all ^-dependence, c/>rc' s (£, C) ~ 
</> r ° es (£) = 0f-(£,O), and obtain 



d i (AN (l d^ cs )--d sW ° + 



1 + V 



dl s AN Cl = 0, 



where 



AiV Cl (0 



«(*) 



exp 



4g(f) 



Integrating Eq. dC8> with the boundary condition 



w 



(C8) 



(C9) 



at 



|£| — > oo and expressing w° in terms of < 
(IB13b> . we obtain 



res 



0«C + 



1 



iraAN^ 



C(x) 



according to 



2(i + VM*) ' 

(C10) 

Renormalizing the potential and coordinate as follows, 



we arrive at an integrodifferential equation 



(Cll) 



d v (p + 



2pe 



V74 



dr 



7I" 3/2 J-oc T -V 

which depends only upon one dimensionless parameter 



P(t) = 



2aAN Cl {0) 



1 



7 



(C12) 



Po 



(C13) 
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Here 



Po = p(to) 



26on (l + V)' 



(C14) 



Multiplying Eq. dC12l by exp(— 77 2 /4), applying the Hilbert 
transform, Eq. (IB 1 3i . and using the identity 



V 



1 



fix) 



dx)dT= -7T 95(77), 



we arrive at a different form of integrodifferential Eq. (IC12I . 



p<p(v) 



1 



2^ 



T — Tj 



e - T2 / 4 d T = S(r}). (C15) 



The function S(r)) in the RHS of Eq. ( IC15I stems from the 
integration 



s(v) = - 



1 



20? 



-r 2 /4 



T — Tj 



dr 



(C16) 



and can be recast as 

S(rj) = rje 



, ,''/ 2 a 

-"/ 4 / e r dr-l 



(C17) 



where erf (a;) denotes the standard error-function. 

Now we discuss some properties of the function S(r)), 
which will be used in Sect. IV Bl and others. Firstly, we observe 
that Eq. JC16i shows that S(rj) a nd ( v /7r/2)?7exp(-?7 2 /4) 
form a Hilbert transform pair, Eqs. JB13h .b). so that 



S(r) ^ 

or = — - — tj exp 



1-coT-r) 2 
Considering the double integral 



>r_ 

4 



(CI 8) 



K = 



dr)V 



30 F(y) 

ooV-V 



with an arbitrary integrable function F(y) and changing the 
order of integration, we verify that K = 0. Applying this to 
Eq. (IC 161 . we conclude that the even function S(rj) obeys 



S{r]) d?y = 0. 



(C19) 



The function S (77) has important power-series approximations 
at sufficiently small and large values of \rj\. The function S(rj) 
can be represented as an infinite Taylor series, 



1 00 



_! rf n (n- 1) 
(2n- 1)! 

^6 

2x6 + 2 x 6 x 10 



n=l 

2 „4 



. (C20) 



which is convergent for all values of 7/. A truncated series with 
a few first terms approximates the function at \rj\ < 2 well. At 
> 4, the function S(rj) can be approximated well by an 
asymptotic series, 



A x (2n- 1)! _ 2 
-J (n — 1)!?7 2 ™ ?7 2 



2x6 2x6x10 



(C21) 

Unlike the Taylor series, this asymptotic series is divergent, so 
that the total number of terms n max should not be too large. 
To conclude this Appendix, we calculate the integral 



1 



2Vi 



V 



d T S(r) 
t — rj 



~^dr, 



(C22) 



needed to obtain the first-order correction of the later-stage 
residual potential in Sect. IV B"T1 Exact integration yields 



J(v) 




(C23) 



The function J(rj) looks qualitatively as S(rf) and to a good 
accuracy can be approximated by a simpler expression 



where the constants / and p are given by 
2 



I = 



0.643, 

1/2 



I + tt/2 



1.546. 



(C24) 

(C25a) 
(C25b) 



We have chosen these constants to provide the best parabolic 
fit between J(rj) and J(ry) at 77 = 0. The biggest difference 
between the exact and approximate expressions for J(rj) is 
near the maxima of J(?7) (|?y| w 2), where it reaches about 
0.016 (less than 10%, see Fig. |8j. In all other locations, the 
functions J (77) and J (77) are much closer to each other. 



APPENDIX D: PROPERTIES OF RESIDUAL POTENTIAL 

In this Appendix, we discuss general properties of solutions 
of Eq. ( 1231 . The solution of ip(x) has the following important 
property, 



ip{x)dx = 0, 



(Dl) 



which follows from general properties of analytic functions 
[see Eq. JB1 H I and symmetry ( I19bt . and is necessary for the 
consistency of Cauchy type integrodifferential Eq. d23t . We 
may rewrite this relation in terms of the x component of the 
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electric field E x — —d x ip. Expressing the residual poten- 
tial as <p(x) = J^ E x dx, substituting this in Eq. iDl\ . and 
changing the order of integration, we obtain 



E x {x)xdx = 0. 



(D2) 



In addition, Eqs. J23h .b) impose restrictions on possible 
asymptotic behavior of the residual potential at large dis- 
tances. For \r)\ 3> 1, we expect the function <p{rj) to behave as 



Pi*)) - > \ 3¥ 



fen 

E 



(D3) 



which corresponds to the expansion of the electric potential 
in multipoles. Note that the lowest-order coefficient C\, in 
accord with Eq. (ICllt . relates to the quadrupole coefficient C 
in Eq. (IB5> as C\ = 2(1 + ip)C. The maximum value of k, 
fcmax, is determined by r\ (the asymptotic series may diverge 
as k — > oo). The electric field proportional to d v <p tends to 
zero as |ry| — > oo. This means that the RHS of Eq. d23ai . 
i.e., the term r], should asymptotically match the integral term 
in the LHS of Eq. (I23ai . The exponentially growing factor 
in front of the principal value integral requires the integral 
to decrease as |»;| — > oo faster than |?7| _ " with positive n. 
Analyzing the asymptotic behavior of the integral, one can 
obtain the following restriction, 




0. 



(D4) 



valid for any positive integer m [in the case of m = 0, we 
obtain Eq. (ID1H . From Eq. (ID4> it follows that the residual 
potential cannot be an exponentially decreasing function of r\ 
as 1 77 1 — * 00, but must have a power-law asymptotic behavior 
(as required by the multipole expansion). Indeed, all coeffi- 
cients C n cannot equal zero because no non-zero ip(x) could 
satisfy Eq. JD4> for all positive integer m. 

The self-similar solution (SSS) for the residual potential, 
defined in Appendix lAl (to the accuracy of an arbitrary con- 
stant), in variables defined by Eq. JC1 II can be written as 
ip(0) + rj 2 /2. Introducing the difference between the actual 
potential <p(rj) and the SSS, 



Tj 

5<p(rj) = <p(ij) - ~ y > 
we rewrite Eqs. d23l > as a Hilbert transform pair: 



d n 5<p(r))e~ 



fxp(r)) 



2p 



20F 



3/2 ' J '-oo T - V 

d T 5(p(r) 



dr. 



* dr, 



n 



(D5) 

(D6a) 
(D6b) 



We expect Sip(rj) w — Xrj 2 as |r7| <C 1 with a positive constant 
A and Sip(rj) w — ^7 2 /2 as \r/\ 3> 1,. Considering the limits of 



small and large values of \t)\, we obtain integral relations 

1 f°° 8 t 5v(t) 



<p{0) 



d(p(rj) 



1 



?7=o 



I'd 



d(r] 2 ) 
lim [vi^rj 2 ] = 



Vnp Jo 
2(1 + 7/0 C x 

1 



T 



4 dr, 



d T 5(p(T) _ £_ 
e 4 



dr, 
(D7) 



d T S(p(T) re 



'dr, 



which describe the parabolic behavior of the residual poten- 
tial near its minimum and express the highest-order multipole 
coefficient in terms of deviation from the self-similar poten- 
tial. Equations ( IDH . JD4I >. JD7> . as well as the exact relations 
obtained in the following Appendix, can be used for checking 
the solution. 



APPENDIX E: DENSITY DISTURBANCE INTEGRALS 

In this Appendix, we derive an important exact property of 
density disturbance integrals (DDIs), which have been applied 
for monitoring the simulation accuracy, see Sect. II VI and IVIII 
We define the two DDIs by integrating parallel and perpendic- 
ular to Bq, 



A7V c (£,i) = / An(£,C,iR, (El a) 

J —OO 

/OO 
An(£,C,*)de, (Elb) 
-00 

where An include both the plasma trail density and distur- 
bances of the background plasma. If we solve Eq. (1151 with 
the self-similar solution (SSS) as the initial condition then the 
DDIs multiplied by t, for all coordinates £ or £, are strictly 
conserved and equal to those for the SSS, even though the 
SSS is, in general, invalid. 

To prove this and obtain the explicit analytical expressions 
for the DDIs, we start from Eqs. dClh .b). Adding them, we 
eliminate ^-derivatives of (f> les , 



d t (iAn) 



d ( (£An) + d c ((An) _ t/>V 2 At7. 

2 1 + 7/; 



1-a 2 



a c [(l + An)5 c ^ res ] = 0. (E2) 



Because both the density disturbance and the potential go to 
zero at infinity, we integrate this equation along the coordinate 
£ from —00 to 00 and obtain 



dt (tAN c ) - d e 



i + tp 



8* AN, 



0. 



(E3) 



Similarly, multiplying Eq. (IClb> by a 2 and adding with (IClat . 
we eliminate ^-derivatives of <^ res . After the integration over 
^ we obtain 



dt (tAN ( ) = d c 



(E4) 
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If we use the SSS as the initial condition at t = to in Eqs. JE3I 
and (|EU then [d t (*AA^)] t=t(j = and [d t (tAN c )] t=to = 0. 
The unique solutions of Eqs. JE3> and JE4t will keep these 
relations in their self-similar form, 



AN^,t) 



2C 




x/tt exp 



for all t. This completes the proof. 



exp 



i/> J A 
4 



(E5a) 
(E5b) 



APPENDIX F: EFFECTIVE INTEGRAL OF DENSITY 
DISTURBANCES 

In this Appendix, we calculate the effective integral over 
near zone, as discussed in Sect. IV D 21 If the trail den- 
sity is the major contribution to the integral over the effec- 
ti ve region then AN Cl w A/V Trail (£), where A7V T r a ii(£) = 

^oMt)(n t /t)exp[-e/Mt)} and q Q = i/>/(l + ip) 
[see Eqs. (IC6> and dC9i in Appendix Icl. If, however, we ex- 
tend the integral over £ to the entire half-axis < ( < oo, 
then we will include the entire background plasma distur- 
bances along this half-axis and obtain AA^ = AiVf/2(£) = 
{rioto/t) exp(— £ 2 /4go), where AN^ is the density dis- 
turbance integral (DDI) over the entire £ axis, Eq. ( IE5at . For 
q(t) ~ 1, the difference between A/Vxraii an d AiV^ is signif- 
icant. The two integrals, however, are both Gaussian functions 
of £ with peaks at £ = 0. Their integrals over the entire £- 
axis, i.e., the 2D integrals of the trail density and the total den- 
sity disturbances over the entire half-space £ > 0, are equal, 
J^AJVWde = fZo AN °°<% = 2n^qo(not /t). The 
true function AN^ (£), taken over the effective region within 
the near zone, < £ < £i, should combine the entire con- 
tribution from the trail density with a part of the background 
density disturbances. It is natural to assume that AN^ (£) 
varies between the two Gaussian functions, A/Vrrail(£) and 
AiVoo (£), and can be approximated by another Gaussian func- 
tion with the peak at £ = 0. This function should have the 
same integral over £, so that it is determined by one parameter 
9, 



q 



1/2 



n t 
t 



exp - — 



(Fl) 



The reason why the integral over £ should be the same is that 
at any £ beyond the trail the corresponding integral of the 
background density disturbances along £, Eq. (IE5b> . is ex- 
ponentially small, as confirmed by our simulations. 

The parameter q(t) is similar to q in iVrrail an d is uniquely 
related to the peak value AiV^ 1 (0). From the derivation of 
Appendixlcl it is clear that the adjusted parameter p is related 
to q by the same relation as p to q, Eq. ( I25> , i.e., p — jqt. 
Thus to find p, we need to estimate AN^ (0), 



1/2 



n t ( np t 



1/2 



(F2) 



The part of the corresponding integral stemming from the trail 
density is determined by q or by p — jqt. To find the ad- 
ditional part stemming from the background density distur- 
bances, we use Eq. J75i . where we should replace p by p be- 
cause these disturbances are determined by the residual poten- 
tial. While the constant value of Sn = 1/1(1 + ip)p is easy 
to integrate, it is not so easy to determine the upper limit of 
integration, i.e., the exact value of £i. 

It is clear, however, that £i should lie within the near zone, 
Eq. @DJ, so that Ci = 13/ a, where a = (1 + ifj) 1/2 O , 
Eq. ( IB2I >. determines the typical scale (~ a -1 ) of the resid- 
ual potential variation along £ and /3 is a numerical factor of 
order unity or less. Adding the two contributions and using 
Eq. JF2I . we obtain 



no 



V tp J 



1/2 



n 



1/2 



(3 



Manipulation with Eqs. 
<A) 3/2 e no(*oPo) 1/2 = (1 
Eq. (IF3> as a quadratic equation for 



V tp j ■ 2(1 + v) 3 / 2 e 0/ 5 ■ 

(F3) 

16) and GHJ yields 2(1 + 

- -0)/7 1 ' 2 , so that we can recast 



where 



P 
P 



13 = 




P , (3 



p 2 



0. 



2/3 [Tt 



i + ip V n P 

Solving Eq. ( IF4I >. we obtain 



p(t) 



P(t) 



l-/3± VI -2/3 



(F4) 



(F5) 



(F6) 



Equation (IF6> contains parameters /3 or /3 which are still un- 
known functions of p and tp. From Eq. JF6I it is clear that 
there exists an upper restriction on them (p(t)/jt — > 1 as 
p -> oo), 



A, 



V5F(1 + V) 



(F7) 



If we start from small values of p, as assumed by our theory, 
then we have to choose the solution with the 'plus' sign in 
front of the square root corresponding to p — > p as j3 — ► 0. 
In this case, the ratio p/p decreases with increasing (3, so that 
p cannot be less than p/A. In our simulations, however, we 
found cases when, at least for some time, p followed Eq. JF6I 
with the 'minus' sign in front of the square root, so that p 
was less than p/A. These special cases (ip = 0.05, Ano < 
50) usually start from sufficiently large values of p which are 
beyond the assumptions of our analytic theory. 

While the strongest discrepancy between p and p takes 
place for p ^> 1, the deviation between the two parameters 
starts developing at p ~ 1. Although the above heuristic 
derivation for p is only valid for /)> 1, we can formally ex- 
tend Eqs. JF4> -l lF6l to p ~ 1 and determine the function (3(p) 
numerically using comparison with simulations, see Sect. lVIll 
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For if) > 0.05, to good accuracy, the numerically found cou- the maximum value of (3 max — 1/2. We extrapolate this to 
pled functions f3(p) and (3(p) can be approximated by smaller ip keeping the same p-dependence as in Eq. jF8> . 



^ R> ^ L ' ^ W M + S + n /S- 5(p) ^< 1+ ^ g t g (F9) 

2.2 + p (l + V)(2.2 + p) V 7rp W~ 4 f 9 9 -U ' W 9r9 9-i-^V n • 1 yj 

(F8) 

Note that for ip = 0.05 at large time, (3 asymptotically reached 



4(2.2 + p) ' 2(2.2 + p) V P 
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